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ABSTRACT 

The Planck full mission eosmic mierowave background (CMB) temperature and E-mode polarization maps are analysed to obtain constraints 
on primordial non-Gaussianity (NG). Using three classes of optimal bispectrum estimators — separable template-fitting (KSW), binned, and 
modal — we obtain consistent values for the primordial local, equilateral, and orthogonal bispectrum amplitudes, quoting as our final result from 
temperature alone f^l^^ = 2.5 + 5.7, f^^^^ = -16 + 70, and f°^^ = -34 + 33 (68 % CL, statistical). Combining temperature and polarization data 
we obtain = 0.8 + 5.0, f^^^^ = -4 + 43, and f^^^ = -26 + 21 (68 % CL, statistical). The results are based on comprehensive cross-validation 
of these estimators on Gaussian and non-Gaussian simulations, are stable across component separation techniques, pass an extensive suite of 
tests, and are consistent with estimators based on measuring the Minkowski functionals of the CMB. The effect of time-domain de-glitching 
systematics on the bispectrum is negligible. In spite of these test outeomes we conservatively label the results including polarization data as 
preliminary, owing to a known mismatch of the noise model in simulations and the data. Beyond estimates of individual shape amplitudes, we 
present model-independent, three-dimensional reconstructions of the Planck CMB bispeetrum and derive constraints on early universe scenarios 
that generate primordial NG, including general single-field models of inflation, axion inflation, initial state modifications, models producing parity- 
violating tensor bispectra, and directionally dependent vector models. We present a wide survey of scale-dependent feature and resonance models, 
accounting for the “look elsewhere” effect in estimating the statistical significance of features. We also look for isocurvature NG, and find no 
signal, but we obtain constraints that improve significantly with the inclusion of polarization. The primordial trispectrum amplitude in the local 
model is eonstrained to be = (-9.0 + 7.7) x 10"^ (68 % CL statistical), and we perform an analysis of trispectrum shapes beyond the local case. 
The global picture that emerges is one of consistency with the premises of the ACDM cosmology, namely that the structure we observe today was 
sourced by adiabatic, passive, Gaussian, and primordial seed perturbations. 

Key words, cosmology: cosmic background radiation - cosmology: observations - cosmology: theory - eosmology: early Universe - eosmology: 
inflation - methods: data analysis 
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1. Introduction 

This paper, one of a set associated with the 2015 release of data 
from the Planck^ mission (Planck Collaboration I 2016), de¬ 
scribes the constraints on primordial non-Cans sianity (NG) ob¬ 
tained using the cosmic microwave background (CMB) maps 
from the full Planck mission, including a first analysis of some 
of the Planck polarization data. 

Primordial NG is one of the most powerful tests of infiation, 
and more generally of high-energy early Universe physics (for 
some reviews, see Bartolo et al. 2004a, Liguori et al. 2010, Chen 
2010b, Komatsu 2010, Yadav & Wandelt 2010). In fact, the sim¬ 
plest models of infiation (characterized by a single scalar field 
slowly rolling along a smooth potential) predict the generation 
of primordial fiuctuations that are almost Gaussian distributed, 
with a tiny deviation from Gaussianity of the order of the slow- 
roll parameters (Acquaviva et al. 2003; Maldacena 2003). The 
2013 Planck results on primordial NG are consistent with such 
a prediction, being compatible with Gaussian primordial fiuc¬ 
tuations: the standard scenario of single-field, slow-roll infia¬ 
tion has survived its most stringent test to date. For example, 
in Planck Collaboration XXIV (2014) we obtained = 2.7 + 

5.8, = -42+75, and/j^f ° = -25 + 39 forthe amplitudes of 

three of the most well-studied shapes of primordial NG. On the 
other hand, it is well known that any deviations from the stan¬ 
dard picture of infiation have the potential to produce distinctive 
NG signatures at a detectable level in the CMB anisotropies.^ 
Therefore, as already shown in Planck Collaboration XXIV 
(2014) (see also Planck Collaboration XXII (2014)) improved 
NG constraints allow severe limits to be placed on various 
classes of inflationary models that extend the simplest paradigm, 
in a way that is strongly complementary to the power-spectrum 
constraints (i.e., scalar spectral index of curvature perturbations 
and tensor-to-scalar amplitude ratio). 

One of the main goals of this paper is to improve NG 
constraints using mainly the angular bispectrum of CMB 
anisotropies, i.e., the harmonic transform of the 3-point angu¬ 
lar correlation function. We also investigate higher-order NG 
correlators like the trispectrum. We follow the same notation 
as Planck Collaboration XXIV (2014). The CMB angular bis¬ 
pectrum is related to the primordial bispectrum 

<0(*i)<I>(*2)3>(*3)> = (27T)W\ki+k2 + k3)B^(ki,k2,k3), (1) 

where the field O, related to the comoving curvature perturba¬ 
tion ^ on super-horizon scales by O = (3/5)f, is such that in the 
matter era, and on super-horizon scales, it reduces to Bardeen’s 
gauge-invariant gravitational potential (Bardeen 1980). The bis¬ 
pectrum B^(ki,k 2 , h) measures the correlation among three per¬ 
turbation modes. If translational and rotational invariance are as¬ 
sumed, it depends only on the magnitude of the three wavevec- 
tors. In general the bispectrum can be written as 

,k2,k2)= fn\^F{ki , ^ 2 , ^ 3 ), ( 2 ) 


' Planck (http://www.esa.int/Planck) is a project of the 
European Space Agency (ESA) with instruments provided by two sci¬ 
entific consortia funded by ESA member states and led by Principal 
Investigators from Prance and Italy, telescope reflectors provided 
through a collaboration between ESA and a scientific consortium led 
and funded by Denmark, and additional contributions from NASA 
(USA). 

^ We refer the reader to Planck Collaboration XXIV (2014) and ref¬ 
erences therein for a detailed summary of the models and underlying 
physical mechanisms generating various types of primordial NG. 


where we have introduced the dimensionless “nonlinearity pa¬ 
rameter” /nl (Gangui et al. 1994; Wang & Kamionkowski 2000; 
Komatsu & Spergel 2001; Babich etal. 2004), measuring the 
NG amplitude. The bispectrum is obtained by sampling triangles 
in Fourier space. The dependence of the function F(ki, ^ 2 , k^) on 
the type of triangle (i.e., the configuration) formed by the three 
wavevectors describes the shape (and the scale dependence) of 
the bispectrum (Babich et al. 2004), which encodes much phys¬ 
ical information. Different NG shapes are linked to distinctive 
physical mechanisms that can generate such NG fingerprints in 
the early Universe. 

In this paper the limits on primordial NG are mainly im¬ 
proved through the use of the full mission data, as well as by 
exploiting the polarization information. 

Planck results on primordial NG also provide a reconstruc¬ 
tion of the full CMB bispectrum through different techniques 
(see Sect. 6.2). This complements (and adds to) the extraction 
of single amplitudes /nl for specific bispectrum shapes. Such 
a reconstruction can point to interesting features in the bispec¬ 
trum signal that go beyond the usual standard scale-invariant 
shapes (such as the well known “local” and “equilateral” con¬ 
figurations). 

As we have seen, the Planck 2013 NG paper 
(Planck Collaboration XXIV 2014) significantly improved 
constraints on the standard primordial NG models with scale- 
invariant local, equilateral or orthogonal shapes. The Planck 
NG paper also included constraints from the modal estima¬ 
tor on a variety of other primordial models, including DBI 
infiation, non-Bunch-Davies models (excited initial states), 
directionally-dependent vector infiation models, warm infiation, 
and scale-dependent feature and resonance models. All scale- 
invariant bispectra were strongly constrained, with the possible 
exception of highly flattened non-Bunch-Davies models. On the 
other hand, the preliminary investigation of primordial oscilla¬ 
tory models seemed to be more promising, in that two specific 
feature models appeared to produce fits of some significance. 
One aim of the present work is to expand the detail and scope 
of investigations of feature and resonant models and to examine 
the significance of these results with a more careful analysis of 
the “look elsewhere” effect, through exploring multi-parameter 
results using large ensembles of Gaussian simulations. Also 
we will thoroughly analyse or re-analyse other primordial NG 
signals that are theoretically well-motivated and those which 
have appeared in the literature since the first data release. These 
include primordial NG arising in the context of infiation models 
where vector fields play a non-negligible role or primordial NG 
generated in the tensor (gravitational waves) perturbations. Each 
of these primordial NG signals carry distinctive signatures that 
may have been imprinted at the inflationary epoch, thus opening 
up a new window into the detailed physics of infiation. 

The paper is organized as follows. In Sect. 2 we briefly 
discuss the primordial NG models that we test in this paper. 
Section 3 summarizes the optimal statistical estimators used 
to constrain the CMB bispectrum and trispectrum from Planck 
temperature and polarization data. In Sect. 4 we discuss the non- 
primordial contributions to the CMB bispectrum and trispec¬ 
trum, including foreground residuals after component separation 
and focusing on the /nl bias induced by the ISW-lensing bispec¬ 
trum. We also analyse the impact on primordial NG estimation 
from the residuals of the deglitching processing. Section 5 de¬ 
scribes an extensive suite of tests performed on realistic simula¬ 
tions to validate the different estimator pipelines, and compare 
their performance. Using simulations, we also quantify the im¬ 
pact on /nl of using a variety of component-separation tech- 
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niques. In Sect. 6 we derive constraints on /nl for the local, 
equilateral, and orthogonal bispectra and present a reconstruc¬ 
tion of the CMB bispectrum. We also present a reconstruction 
of the primordial curvature fluctuations. In Sect. 7 we validate 
these results by performing a series of null tests on the data to 
assess the robustness of our results. Section 8 investigates scale- 
dependent NG models and other selected bispectrum shapes. 
Section 9 presents the Planck limits on the CMB trispectrum. 
In Sect. 10 we provide constraints on CMB local bispectrum 
and trispectrum from Minkowski functionals. In Sect. 11 we dis¬ 
cuss the main implications of Planck's constraints on primor¬ 
dial NG for early Universe models. We conclude in Sect. 12. 
Appendix A contains the derivation of the statistical estima¬ 
tor for the amplitudes characterizing a “direction-dependent” 
primordial non-Canssianity. Appendix B contains some details 
about Minkowski functionals. 


2. Models 

In this Section we briefly highlight the classes of inflationary 
models investigated in this paper, and describe the distinctive 
NG they generate. Within each class a common underlying phys¬ 
ical process gives rise to the corresponding NG shape, illustrated 
by concrete realizations of inflationary models. For each class 
we therefore provide the explicit form of the bispectrum shapes 
chosen for the data analysis, emphasizing extensions with vari¬ 
ants and distinctly new shapes beyond those already described 
in Planck Collaboration XXIV (2014). 


2.1. General single-field models of inflation 


This class of models includes inflationary models with a non¬ 
standard kinetic term (or more general higher-derivative inter¬ 
actions), in which the inflaton fluctuations propagate with an 
effective sound speed Cs which can be smaller than the speed 
of light. For example, models with a non-standard kinetic term 
are described by an inflaton Lagrangian X = P(X, 0 ), where 
X = g^^d^(p dy(p, with at most one derivative on 0, and the sound 
speed is given by c^ = {dPjdX) j{dPjdX + 2X{d^PldX^)). 

The NG parameter space of this class of models is generi- 
cally well described by two NG shapes — “equilateral” and “or¬ 
thogonal” (Senatore et al. 2010) — since usually there are two 
dominant interaction terms of the inflaton field giving rise to the 
overall NG signal. One of these typically produces a bispectrum 
very close to the equilateral type with /nl ~ in the limit 
Cs 1 (Chen et al. 2007b; Senatore et al. 2010). 

The equilateral-type NG is well approximated by the tem¬ 
plate (Creminelli et al. 2006) 




T 4-ns 

/C 2 


T 4-ns 
^2 ^3 


T 4-ns 

/C^ /Cj 


(kik2k3)^^^ "0/3 


1 


jM-ns)/3j2(4-n,)/3j4-n, 


-r (5 permutations) 


}■ 


( 3 ) 


Here P^(k) = A/k^~^^ represents Bardeen’s gravitational po¬ 
tential power spectrum, being the normalization and n^ the 
scalar spectral index. DBI inflationary models based on string 
theory (Silverstein & Tong 2004; Alishahiha et al. 2004) pro¬ 
vide physically well-motivated examples of the P(X, 0)-model. 
They are characterized by an almost equilateral NG with “ 
-(35/108)c“^ for Cs 1, which typically is < -5. 


The “orthogonal” shape template is (Senatore et al. 2010) 


BZ^(kuk2,k,) = 6A^fX' 
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( 4 ) 


Equilateral and orthogonal shapes emerge also in models char¬ 
acterized by more general higher-derivative interactions, such 
as ghost inflation (Arkani-Hamed et al. 2004), effective field 
theories of inflation (Cheung et al. 2008; Senatore et al. 2010; 
Bartolo et al. 2010a), or the so “Galileon-like” models of infla¬ 
tion (see, e.g., Burrage et al. 2011). The latter model is con¬ 
structed starting from some specific underlying symmetry for 
the inflaton field, and is characterized by strongly constrained 
derivative interactions. 


2.2. Multi-field models 


This class of models is characterized by the presence of addi¬ 
tional light scalar degrees of freedom besides the inflaton, whose 
fluctuations give rise, or contribute, to the final primordial cur¬ 
vature perturbation at the end of inflation. This includes the case 
of “multiple-field inflation”, where inflation is driven by more 
than one scalar held, as well as scenarios in which additional 
scalar fields remain subdominant during the inflationary expan¬ 
sion. From the point of view of primordial NG, the element in 
common to all these models is that a potentially detectable level 
of NG in the curvature perturbation is generated via a transfer 
of super-horizon non-Gaussian isocurvature perturbations in the 
second held (not necessarily the inflaton) to the adiabatic (cur¬ 
vature) density perturbations, accompanied by nonlinearities in 
the transfer mechanism. This process typically takes place on 
super-horizon scales, thus implying a local form of NG in real 
space. When going to Fourier space, this leads to a correlation 
between large and small scale modes. The bispectrum for this 
class of models is indeed largest on so-called “squeezed” trian¬ 
gles (^1 k 2 - ^ 3 ). The local bispectrum is (Falk et al. 1993; 
Gangui et al. 1994; Gangui & Martin 2000; Verde et al. 2000; 
Wang & Kamionkowski 2000; Komatsu & Spergel 2001) 




•local 


ihMM) = 2f^C'[P^{h)P^{k2) + P^{h)P^(k3) 




0^2 /-local 

/nl 


1 


T 4-ns i4 

IK 1 IKr% 


cycl. 


( 5 ) 


There is a broad literature on examples and specific realiza¬ 
tions of this transfer mechanism from isocurvature to adi¬ 
abatic perturbations (Bartolo et al. 2002; Bernardeau & Uzan 
2002; Vemizzi & Wands 2006; Rigopoulos et al. 2006, 2007; 
Lyth & Rodriguez 2005; Tzavara & van Tent 2011; for a review 
on NG from multiple-field inflation models, see Byrnes & Choi 
2010). An alternative, important possibility is the cur- 
vaton model (Mollerach 1990; Linde & Mukhanov 1997; 
Enqvist & Sloth 2002; Lyth & Wands 2002; Moroi & Takahashi 
2001). In this type of scenario, a second light scalar held, sub¬ 
dominant during inflation, decays after inflation, generating pri¬ 
mordial density perturbations with a potentially high level of 
NG (e.g., Lyth & Wands 2002; Lyth et al. 2003; Bartolo et al. 
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2004c). In the (simplest) adiabatic curvaton models, the local 
/nl parameter was found to be (Bartolo et al. 2004c, b) = 

(5/4rD)-5rD/6-5/3, when the curvaton field has a quadratic po¬ 
tential (Lyth & Wands 2002; Lyth et al. 2003; Lyth & Rodriguez 
2005; Malik & Lyth 2006; Sasaki et al. 2006). In the previous 
formula, Td = [ 3 pcurvaton/( 3 pcurvaton + 4 pradiation)]D IS the CUrva- 
ton decay fraction” evaluated at the epoch of the curvaton decay 
in the sudden decay approximation. It is then easy to see that, for 
low values of td, a high level of NG can be generated. ^ 

2.3. Isocurvature non-Gaussianity 

Isocurvature NG, which was only sketched from the purely the¬ 
oretical point of view in the 2013 paper, can now be analysed 
thanks to the polarization information. 

In most of the models mentioned above, the main focus is on 
the level of primordial NG in the final curvature perturbation f. 
However, in inflationary scenarios where different scalar fields 
play a non-negligible role, residual isocurvature perturbation 
modes can remain after inflation. Isocurvature modes are usu¬ 
ally investigated by considering their contribution to the power 
spectrum. However, if present, they would also contribute to the 
bispectrum, producing in general both a pure isocurvature bis¬ 
pectrum and mixed bispectra because of the cross-correlation 
between isocurvature and adiabatic perturbations (Komatsu 
2002; Bartolo et al. 2002; Komatsu et al. 2005; Kawasaki et al. 
2008; Langlois et al. 2008; Kawasaki et al. 2009; Hikage et al. 
2009 ; Langlois & Lepidi 2011; Langlois & van Tent 2011; 
Kawakami et al. 2012; Langlois & van Tent 2012). While one 
might expect isocurvature NG to be negligible, since both (lin¬ 
ear) isocurvature modes and (adiabatic) NG appear to be very 
small, and searches for isocurvature NG using WMAP data did 
not lead to any detections (Hikage et al. 2013a,b), this expec¬ 
tation can be tested at significantly higher precision by Planck. 
Moreover, there exist inflation models (Langlois & Lepidi 2011) 
where isocurvature modes, while remaining a small fraction in 
the power spectrum, would dominate the bispectrum. 

At the time of recombination there are in principle four pos¬ 
sible distinct isocurvature modes (in addition to the adiabatic 
mode): cold dark matter (CDM); baryon; neutrino density; and 
neutrino velocity isocurvature modes (Bucher et al. 2000). In 
this paper we will only consider isocurvature NG of the lo¬ 
cal type and always limit ourselves to considering the adia¬ 
batic mode together with just one type of isocurvature mode 
(considering each of the four types separately). Otherwise the 
number of free parameters becomes so large that no mean¬ 
ingful limits can be derived. Moreover, we assume the same 
spectral index for the primordial isocurvature power spectrum 
and the isocurvature-adiabatic cross-power spectrum as for the 
adiabatic power spectrum. Under those assumptions, as shown 
by Langlois & van Tent (2011), we have in principle six inde¬ 
pendent /nl parameters: the usual purely adiabatic one; a purely 
isocurvature one; and four correlated ones. 


^ NG perturbations can arise also at the end of inflation, e.g., 
from nonlinearities during the (p)reheating phase (e.g., Enqvist et al. 
2005; Chambers & Rajantie 2008; Barnaby & Cline 2006; see 
also Bond et al. 2009) or from fluctuations in the inflaton decay rate 
or interactions, as found in modulated (p)reheating and modulated 
hybrid inflation (Kofman 2003; Dvali et al. 2004a,b; Bemardeau et al. 
2004; Zaldarriaga 2004; Lyth 2005; Salem 2005; Lyth & Riotto 2006; 
Kolb et al. 2006; Cicoli et al. 2012). 


The primordial shape templates are a generalization of 
Eq. (5), see Langlois & van Tent (2011, 2012): 

B"’^{kuk2,h) = 2fl^l’^P,t,{k2)P^{h) + 2f^‘P^{h)P,t,{k3) 

+2f^l-’PMP<t>ik2), (6) 

where /, /, K label the different modes (adiabatic and isocurva¬ 
ture). The invariance under the simultaneous exchange of two 
of these indices and the corresponding momenta means that 
hence reducing the number of independent param¬ 
eters from eight to six, in the case of two modes. The different 
bispectra vary most importantly through the fact that different 
types of radiation transfer functions g^^{k) are used to project the 
primordial template onto the CMB: the reduced bispectra are of 
the form 



/^CO 

(7) 

with 



a‘f{r) = 

Cdkje{kr)g\{k), 

(8) 

P'eir) = 

^ J'Cdkjt{kr)g[{k)Pq>{k). 

(9) 


Here is the spherical bessel function and we use the notation 
(^i^ 2 ^ 3 ) = [G^ 2^3 + 5perm.]/3!. In addition to the isocurvature 
index, each transfer function carries a polarization index that we 
do not show here. It is important to note that, unlike the case of 
the purely adiabatic mode, the inclusion of polarization improves 
the constraints on the isocurvature NG significantly, as predicted 
by Langlois & van Tent (2011, 2012). 

2.4. Resonance and axion monodromy models 

Oscillatory models for NG are physically well-motivated. Large- 
field inflation faces an inherent UV completion problem because 
the inflaton field is required to move over large distances in field 
space relative to the Planck mass mpi. An effective shift sym¬ 
metry can enforce potential flatness and this can be naturally 
implemented in a string theory context with axions and a period¬ 
ically modulated potential, so-called “axion monodromy” mod¬ 
els. This periodicity can generate resonances in the inflationary 
fluctuations with logarithmically-spaced oscillations, creating 
imprints in the power spectrum, the bispectrum and trispectrum 
(Chenetal. 2008; Flaugeretal. 2010; Hannestad et al. 2010; 
Flauger & Pajer 2011). On the other hand, sharp features or 
comers in an inflationary potential can temporarily drive the 
inflaton away from slow-roll; these large changes in the field 
and derivatives can create evenly-spaced oscillations, to be dis¬ 
cussed in the next subsection. However, in multifield models 
residual oscillations after corner-turning can also lead to log¬ 
spaced oscillations, just as in the resonance models (Chen 2011; 
Achucarro et al. 2011; Battefeld et al. 2013; Chenetal. 2015). 
A preliminary search for bispectmm resonance signals was per¬ 
formed in the first Planck analysis (Planck Collaboration XXIV 
2014) and our purpose here is to substantially increase the fre¬ 
quency range and number of models investigated. 

Simple resonance model: Periodic features in the inflationary po¬ 
tential can induce oscillations with frequency oj that can resonate 
through any interactions with the inflationary fluctuations, con¬ 
tributing to the bispectrum. Provided that oj > H, this mode 
starts inside the horizon but its frequency decreases as it is 
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stretched by inflation, until frozen when oj ^ H. Thus peri¬ 
odic features introduce a driving force which can scan across 
a wide range of frequencies. The simplest basic behaviour of 
such resonant models yields logarithmic stretching and can be 
described by the non-scale-invariant shape (see Chen et al. 2008; 
Chen 2010b) 


^^2 ^es 

B%\ki, k2, h) = +k2 + h) + (p] , (10) 

{kik2k2Y 

where the constant C = l/ln(3^:,), K is a wavenumber associ¬ 
ated with the periodicity, and 0 is a phase. These oscillations 
constructively and destructively interfere with the oscillations 
created by the CMB transfer functions, introducing additional 
nodal points in the CMB bispectrum. 

Generalized resonance models: In a more general context, it 
is possible to have more complicated resonant shapes and en¬ 
velopes. Resonant single-fleld models with varying sound speed 
Cs generate three leading-order bispectrum terms (Chen 2010a): 

B^^^-^\kuk2,h) = sin[C \n(ki +k2 + h) + cf] 

+ cos {C\n{h +k 2 + kY + <!>] 

+ sin [C \niki +k2 + k^ + cf>]}. (11) 

The first term on the right-hand side of Eq. (11) is the basic res¬ 
onant shape given in Eq. (10), while the second and third terms 
have the same oscillatory behaviour, but modulated by a (mildly) 
flattened shape, and an equilateral shape respectively. The third 
term is in fact the second generic shape arising in effective held 
theory and correlates well with the equilateral shape in Eq. (12). 
The second term in Eq. (11) weakly favours flattened triangles, 
but there are regimes for resonant models that can generate much 
stronger flat shapes. If the resonance begins very deep inside the 
horizon, then the second (negative energy) mode can also make 
a significant contribution that is associated with enfolded or flat 
bispectra; this is similar to having an excited initial state or non- 
Bunch-Davies (NBD) vacuum. 

With these two physical motivations in mind we also inves¬ 
tigate classes of models with resonant oscillations modulated by 
both the equilateral and flattened shapes, defined by 

S^Hkuk2,k3)=^^, = (12) 


where ki = ^ 2+^3 (here, for simplicity we ignore the spectral 
index dependence of the equilateral shape in Eq. 3). The corre¬ 
sponding equilateral and flattened resonant bispectra ansatze are 
then 


^res = S^^(ki,k 2 , ^3) X ^3) 

= ,, {’^^32 FFF +k2 + k3) + (p] , ( 13 ) 

\k 1 k 2 k 23 ) k\k 2 k^ 

^res-flat(^^^ ^ 2 , fcj) s ^ 2 , ^ 3 ) X B’^^\kuk2, h) . (14) 

We note that typically non-Bunch-Davies bispectra can be 
much more sharply peaked in the flattened or squeezed limits 
than Eq. (14), but our purpose here is to determine if this type of 
resonant model is favoured by the Planck data; that is, whether 
Eq. (14) warrants further investigation with other flattened pro¬ 
files. 


2 . 5 . Scale-dependent oscillatory feature models 

Temporary violations of slow-roll inflation can occur if there are 
sharp features in the inflationary potential (Chen et al. 2007a), 
as well as changes in the sound speed Cs or sharp turns in held 
space in multifield inflation. The inflaton held makes temporary 
departures from the attractor solution, which typically have a 
strong scale-dependent running modulated by a sinusoidal os¬ 
cillation; there are model-dependent counterparts in the power 
spectrum, bispectrum, and trispectrum. Eor example, sharper or 
narrower features induce a relatively larger signal in the bispec¬ 
trum (see e.g., Chen 2010b). An example is the analytic envelope 
solutions predicted for both the power spectrum and bispectrum 
for the single held models with a specific inflaton feature shape 
(Adshead et al. 2012); a search for these was presented previ¬ 
ously in the Planck Inflation paper (Planck Collaboration XXII 
2014) and likewise no significant signal was found using the cor¬ 
responding bispectrum envelopes at the available modal reso¬ 
lution (Planck Collaboration XXIV 2014). In this new analysis, 
we will emphasize the search for generic oscillatory behaviour 
in the data over a larger range in modal resolution, although we 
will also look for the shapes predicted for simple features in sin¬ 
gle held models. 

Constant feature model: In the previous investigation of Planck 
data using a coarse parameter grid (Planck Collaboration XXIV 
2014), we searched for the simplest ansatz for an oscillatory bis¬ 
pectrum signal (Chen et al. 2007a): 

^^2 ^feat 

, k2, ks) = sin [w(fci + ^2 + ^3) + ^] , ( 15 ) 

{kik2k3y 


where 0 is a phase factor and tu is a frequency associated with 
the specific shape of the feature in the potential that disrupts the 
slow-roll evolution. In the earlier analysis, we also considered 
a damping envelope, which slightly increased the apparent sig¬ 
nificance of the best-fit feature models, though at the cost of an 
additional parameter (see single-fleld solutions below). 

Generalized feature models: Here, we again search for oscil¬ 
latory signals in a model-independent manner. We will modu¬ 
late the bispectrum cross-sections with the physically motivated 
equilateral and flattened shapes, reflecting the physical contexts 
in which they could have been generated, as for the resonant 
models described above in Eq. (11). If there are potential fea¬ 
tures in a model with a varying sound speed, then we can expect 
there to be oscillatory contributions to the bispectrum signal with 
a dominant equilateral shape. Motivated by the equilateral reso¬ 
nance model in Eq. (11), we will search for the following equi¬ 
lateral feature ansatz: 


^feat ^k 2 ,kf} = S ®^(^i, ^2, ^3) X , ^2, ^3) 


6A^C-^hk2h 

(^i^2^3)^ hhh 


sin \ojiki -I- ^2 + ^ 3 ) + 0] • 


(16) 

(17) 


Eor extremely sharp features, it is possible to excite the inflation¬ 
ary fluctuations as if there were a non-Bunch Davies vacuum: 
the oscillatory signal becomes modulated with a flattened shape 
(Chen et al. 2007a). Again, motivated by the enfolded resonance 
model in Eq. (14), we take the following simple flattened ansatz: 


^feat , ^ 2 , ^ 3 ) X , ^ 2 , ^ 3 ) • (18) 


Although the exact profile of the flattened shape can be much 
more highly peaked on the faces in these NBD models, this 
ansatz should be adequate for testing whether these models are 
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favoured. We note that while the power spectrum is insensi¬ 
tive to the underlying scenario creating the features, the bispec¬ 
trum shape will reveal whether features arise from varying sound 
speed or highly excited features in the potential. 

Single field feature solutions: Here we use the full analytic bis¬ 
pectrum solution given by Adshead et al. (2012), but the domi¬ 
nant leading-order behaviour takes the form 

^^2 cos 

, k2, h) = —^^K^DiacoK) cos{ojK) , (19) 

{k\k2hy 

where K = k\+k 2 +k 2 > and DiacoK) = aojfiK sinh(Qf 6 t;^)) is an 
envelope function, with parameter a setting an overall cut-off for 
the bispectrum at large wavenumbers or multipoles. This enve¬ 
lope and the overall scaling distinguishes this realistic case 
from the simple separable constant feature ansatz of Eq. (15). 
We shall allow the envelope parameter a to vary from a = 0, 
with no envelope (the infinitely thin limit for a feature in the 
potential) through to large a, with a narrow domain for the bis¬ 
pectrum. Alternative analytic solutions where the bispectrum is 
created by a variation in the sound speed are dominated by the 
K sin(tu^) term, as in 

^^2 sin 

B^^^{kuk2M) = —^^KD{acjK)^,m{cjK ). ( 20 ) 

{kik2k2y 

For the simplest models there is a predicted relationship be¬ 
tween the power spectrum and bispectrum amplitude (e.g., see 
also Achucarro et al. 2013 for a two-field model). We note that 
typically the power spectrum has larger signal-to-noise at low 
frequency (i.e., below oj ^ 1000 ) while the bispectrum domi¬ 
nates at higher frequency. 


2.6. Non-Gaussianity from excited initial states 


It is well known that if the initial vacuum state for infla¬ 
tion is excited and deviates from the standard Bunch-Davies 
vacuum, then measurable non-Gaussianities can be produced 
(Chenetal. 2007b; Holman & Tolley 2008; Meerburg et al. 
2009; Ashoorioon & Shiu 2011). These models generically lead 
to non-Gaussianity that peaks in the flattened limit, where ki -\- 
k 2 ~ ^ 3 , and also often has oscillatory behaviour. Here we 
constrain the same selection of templates found in the 2013 
Planck analysis, namely the fiat model in Eq. (12), Non-Bunch- 
Davies (NBD) (Chenetal. 2007b), NBDl and NBD2 models 
(Agullo & Parker 2011) (now called “NBDl cos” and “NBD2 
cos”) and NBD3 (Chen 2010b). We also introduce three new 
templates, NBD sin which is motivated by (Chen 2010a) and 
takes the form 


fiNED-sin (;ti,jt2,A:3) 


9^2 yi-NED-sin 

-^NL ( -coh , „-0jh 

{kik2k-if 




X sin {(jjK -I- 0) , 


( 21 ) 

where again ^ -r ^2 + ^3 and ki = K - 2ki. The other two 
templates are extensions of the NBDl cos and NBD2 cos models 
found in Agullo & Parker (2011) and take the form 

2A^ ^NBD/-sin 

{kik2hY 

X sin( 6 c)^i)/^i -r 2perm.], (22) 

where f\ik\\k 2 ,ky) = k\{k^ k'^)l2, which is dominated by 

squeezed configurations, and / 2 (/^i;^ 2 ,^ 3 ) = ^ 2 ^ 3 , which has a 
flattened shape. 


2.7. Directional-dependence motivated by gauge fields 

Some models where primordial vector fields are present during 
inflation predict interesting NG signatures. This is the case of 
a coupling of the infiaton field ip to the kinetic term of a gauge 
field A^, X contains where = d^Ay - dyA^ and 

the coupling Y(p)F^ is chosen so that scale invariant vector per¬ 
turbations are produced on superhorizon scales (Bamaby et al. 
2012b; Bartolo et al. 2013a). The bispectrum turns out to be the 
sum of two contributions: one of the local shape; and another that 
is also enhanced in the squeezed limit (^1 ^2 - ^ 3 ), but fea¬ 

turing a non-trivial dependence on the angle between the small 
and the large wave vectors through the parameter pu = ki • ^2 
(where k = kjk) as /i^ 2 - Also, primordial magnetic fields sourc¬ 
ing curvature perturbations can cause a dependence on both p 
and p^ (Shiraishi et al. 2012). 

We can parametrize these shapes as variations on the local 
shape (Shiraishi et al. 2013a), as 

Bis>(kuk 2 ,k 2 ) = ^CL[PL(jJi 2 )Pis>(ki)P^(k 2 ) + 2perm.], (23) 

L 


where Pl{p) is the Legendre polynomial with P{) = Pi = p, 
and P 2 = \Op^ - 1). For example, for L = 1 we have the shape 


'^{ki,k2,h) = 




(kik2k3y 


-^{k\ +kl- kj) + 2 perm. 


.(24) 


The local template corresponds to c/ = 2 /nl^/o. Here and in the 
following the nonlinearity parameters are related to the cl 
coefficients by cq = 2 /^°, ci = -4 /^l'’ ^nd C 2 = -16/^^. 
The L = 1,2 shapes exhibit sharp variations in the flattened 
limit, for example for -l- ^2 ~ ^ 3 , while in the squeezed limit, 
L = 1 is suppressed whereas L = 2 grows like the local bispec¬ 
trum shape (i.e., the L = 0 case). The fi(p)F^ models predict 
C 2 = Co 12, while primordial curvature perturbations sourced by 
large-scale magnetic fields generate non-vanishing cq, ci, and 
C 2 . Quite interestingly, in the proposed “solid inflation” scenario 
(Endlich et al. 2013; see also Bartolo et al. 2013b; Endlich et al. 
2014; Sitwell & Sigurdson 2014; Bartolo et al. 2014) bispectra 
similar to Eq. (23) can be generated, in this case with C 2 » 
Co (Endlich et al. 2013, 2014). Therefore, measurements of the c/ 
coefficients can be an efficient probe of some detailed aspects of 
the inflationary mechanism, such as the existence of primordial 
vector fields during inflation (or a non-trivial symmetry structure 
of the infiaton fields, as in solid inflation). 


2.8. Non-Gaussianity from gauge-field production during 
axion inflation 

The same shift symmetry that leads to axion (monodromy) 
models of inflation (Sect. 2.4) naturally allows (from an ef¬ 
fective field theory point of view) for a coupling between a 
pseudoscalar axion infiaton field and a gauge field of the type 
X ^ -{alAf)(pF^^F^y, where the parameter a is dimensionless 
and / is the axion decay constant {F^^ = This sce¬ 

nario has a rich and interesting phenomenology both for scalar 
and tensor primordial fluctuations (see, e.g., Barnaby & Peloso 
2011; Sorbo 2011; Bamaby et al. 2011, 2012c; Linde et al. 
2013; Meerburg & Pajer 2013; Ferreira & Sloth 2014). Gauge 
field quanta are produced by the background motion of the in¬ 
fiaton field, and these in turn source curvature perturbations 
through an inverse decay process of the gauge field. A bispec- 
tmm of curvature fluctuations is generated as (Bamaby et al. 
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2011; Meerburg & Pajer 2013)'^ 


^inv.dec 


= 6A^ 


^inv.dec 

^NL 


Wk] m.AA) ’ 


(25) 


(5/3)^= (5/3)^5^^^V/nl^^^) in the equilateral limit, yield¬ 
ing 


/•tens 

^NL 


5+++(A:i,A:2,A:3) 

lim-n-• 

k,^k F^''\kuk2,h) 


(28) 


where the exact expression for the function /s can be 
found in equation (3.29) of Barnaby et al. (2011) (see 
also Meerburg & Pajer 2013). Here ^ characterizes the coupling 
strength of the axion to the gauge field ^ = a\^\l{2fH). The in¬ 
verse decay bispectrum peaks for equilateral configuration, since 
dip is mostly sourced by the inverse decay (6A -\-6A ^ 6p), when 
two modes of the vector fields are of comparable magnitude (the 
correlation with the equilateral template is 94 % and with the or¬ 
thogonal one is 4 %). We do however constrain the exact shape in 
Eq. (25), without resorting to the equilateral template. Another 
interesting observational signature that can shed light on the role 
played by pseudo-scalars in the early Universe is provided by 
tensor NG, to which we turn next. 

2.9. Parity-violating tensor non-Gaussianity motivated by 
pseudo-scalars 

While the majority of the studies on primordial and CMB 
NG focus on the scalar mode, tensor-mode NG has been at¬ 
tracting attention as a probe of high-energy theories of grav¬ 
ity (e.g., Maldacena & Pimentel 2011; McFadden & Skenderis 
2011; Sodaetal. 2011; Shiraishi et al. 2011; Gaoetal. 2011) 
or primordial magnetic fields (Brown & Crittenden 2005 ; 
Shiraishi et al. 2012; Shiraishi 2012).^ 

Recently, the possibility of observable tensor bispectra has 
been vigorously discussed in a model where the infiaton couples 
to a pseudoscalar field (Barnaby et al. 2012a; Cook & Sorbo 
2013; Ferreira & Sloth 2014). In this model, through the grav¬ 
itational coupling to the U(l) gauge field, gravitational waves 
(hij = I = 2^=+receive NG corrections, where 
only one of the two spin states is enhanced. The bispectrum, 
generally formed as 

|l~| = i2nfd^^\ki +k2 + k3)Bl''^'\ku k2, ks) ,(26) 

is accordingly polarized, with » This 

NG enhancement is a sub-horizon effect and therefore is 
maximized at the equilateral limit (ki ^ ^2 - ^ 3 ) (Cook & Sorbo 
2013). 

A model-independent template of the equilateral-type polar¬ 
ized tensor bispectrum is given by (Shiraishi et al. 2013b, 2015) 

, (27) 

with the polarization tensor obeying \k) = 2ds-s' 

and = e^y\k) = e^^^{-k). We here have introduced 

a tensor nonlinearity parameter, by normalizing with the equi¬ 
lateral bispectrum template of curvature perturbations = 

^ For simplicity we assume a scale-invariant bispectrum. 

^ See Planck Collaboration XIX (2016) for the Planck constraints on 
magnetically-induced NG. 


The template Eq. (27) can adequately reconstruct the ten¬ 
sor bispectra created in the pseudoscalar inflation models^ 
(Shiraishi et al. 2013b), and thus the amplitude is directly 
connected with the model parameters, e.g., the coupling strength 
of the pseudoscalar field to the gauge field ^ (for details see 
Sect. 11). 

The CMB temperature and £’-mode bispectra sourced by the 
parity-violating tensor NG have not only the usual parity-even 
(G+^ 2+^3 = even) signals but also parity-odd (G+^ 2+^3 = odd) 
contributions, which cannot be sourced by known scalar bispec¬ 
tra (Kamionkowski & Souradeep 2011; Shiraishi et al. 2011). 
Moreover, their shapes are mostly distinct from the scalar tem¬ 
plates, due the different radiation transfer functions; hence they 
can be measured essentially independently of the scalar NG 
(Shiraishi et al. 2013b). The analysis of the WMAP temperature 
data distributed in -r G + G = odd configurations leads to an 
observational limit= (0.8+1.1)xl0"^ (Shiraishi et al. 2015). 
This paper updates the limit, by analysing both parity-even and 
parity-odd signals in the Planck temperature and £’-mode polar¬ 
ization data. 

3. Statistical estimation of the CMB bispectrum for 
polarized maps 

We now provide a brief overview of the main statistical tech¬ 
niques that we use to estimate the nonlinearity parameter /nl 
from temperature and polarization CMB data, followed by a de¬ 
scription of the data set that will be used in our analysis. 

The CMB temperature and polarization fields are character¬ 
ized using the multipoles of a spherical harmonic decomposition 
of the CMB maps: 

AT j 

— (h) = 2_^af^Ye„,ih), 

£m 

E(h) . (29) 

£m 

At linear order, the relation between the primordial perturbation 
field and the CMB multipoles is (e.g., Komatsu & Spergel 2001) 

r (Pk 

«i=4^(-0G -^<^(k)Y(m(m^e(k), (30) 

where X = {T,E} denotes either temperature or T-mode polar¬ 
ization, O is the primordial gravitational potential, and Aj rep¬ 
resents the linear CMB radiation transfer function. 

The CMB angular bispectrum is the three-point correlator of 
the aftn^ 

(31) 

where Xf = {T,E}. If the CMB sky is rotationally invariant, and 
the bispectra we are considering have even parity (which is true 
for combinations of T and E), then the angular bispectrum can 
be factorized as 

^ The form of the tensor bispectrum is the same whether the infiaton 
field is identified with the pseudoscalar field or not. 
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where is the so-called reduced bispectrum, and 

is the Gaunt integral, defined as the integral over the solid angle 

of the product of three spherical harmonics, 

f (33) 

The Gaunt integral (often written in terms of Wigner ?>j- 
symbols) enforces rotational symmetry, and restricts attention to 
a tetrahedral domain of multipole triplets {^ 1 ,^ 2 , ^ 3 }, satisfying 
both a triangle condition and a limit given by some maximum 
resolution ^max (the latter being defined by the finite angular res¬ 
olution of the experiment under study). 

Our goal is to extract the nonlinearity parameter /nl from 
the data, for different primordial shapes. To achieve this, we es¬ 
sentially fit a theoretical CMB bispectrum ansatz to the 

observed 3-point function. Theoretical predictions for CMB an¬ 
gular bispectra arising from early Universe primordial models 
can be obtained by applying Eq. (30) to the primordial bispectra 
of Sect. 2, (see e.g., Komatsu & Spergel 2001). Optimized cubic 
bispectrum estimators were introduced by Heavens (1998), and 
it has been shown that for small NG the general optimal polar¬ 
ized /nl estimator can be written as (Creminelli et al. 2006) 


f y y y 
“ A/^ Zj Zj Zj 


XiX ti,mi 


/p-1 ^^3^3/ 7^3 

{31713Am' ' 


(c; 


-1 




, (34) 


where is a suitable normalization chosen to produce unit re¬ 
sponse to • Note that we are implicitly defining a suitable 

normalization convention so that - and 

is the value of the theoretical template when /nl = 1. C“Ms the 
inverse of the block matrix: 





TT 

~E7r 



(35) 


and the blocks represent the full TT, TE, and EE covariance ma¬ 
trices, with being the transpose of All quantities in the 
previous equation (i.e., CMB multipoles, bispectrum template 
and covariances matrices) are assumed to properly incorporate 
instrumental beam and noise. 

As standard for these estimators, we note in square brack¬ 
ets (below) the presence of two contributions. One is cubic 
in the observed a^m^, and correlates the bispectrum of the 
data to the theoretical fitting template • This is generally 
called the “cubic term” of the estimator. The other contribu¬ 
tion is linear in the observed (“linear term”). This part 
corrects for mean-field contributions to the error bars, intro¬ 
duced by rotational invariance-breaking features, such as a mask 
or anisotropic/correlated instrumental noise (Creminelli et al. 
2006; Yadav et al. 2008/ 

The inverse covariance filtering operation implied by 
Eq. (34) is a challenging numerical task, which has been 
successfully performed only recently (Smith et al. 2009; 
Eisner & Wandelt 2012). This step can be avoided by working 
in the “diagonal covariance approximation”. In this approach, 
the estimator is built by neglecting off-diagonal entries of the 
covariance matrix in the cubic term in Eq. (34), and then finding 


the linear term that minimizes the variance for this specific cubic 
statistic. Applying such a procedure yields (Yadav et al. 2007) 




{\ h h /p-1 AAi /p- 


1 A2^2 ("0“ 1 ^2^3 5 th 

T2 ^ T3 ^ 1 ^ 2 ^ 


Xi,X\ {i,mi 


^4^2 ^^3^3 


_ X' _ x;x' x' 

i\m\,i2^2^t3m3 {\m\,{3m3^{2^2 


X'X' 

^2 ^2? A ^3 


a 


{\m\ 


(36) 


where ^ is the inverse of the 2x2 matrix 


This expression can also be written as 


(37) 


(38) 


where the observed (reduced) bispectrum includes the linear cor¬ 
rection term and the inner product is defined as 


{b^,b^) = 




X1X2X3, 

A^ 2 A 


A 


h 


2 

{\t 2 t 2 , 


(C-‘)U‘ 


(39) 

1 X 2 X' 1 X 3 X' x;x'x',fi 

(G )^^ (G )^^ 


XiXi ti 


with 


, 1(2(1 + 1 )( 2/2 + 1)(2/3 + 1 ) ki (2 ( 3 ] 

where the last term is the Wigner 3j symbol. The denominator in 
Eq. (38), (Z?*, Z?*) is the normalization constant N. 

The price to pay for the simplification obtained in Eq. (36) is, 
in principle, loss of optimality. However, in practice we found in 
our previous temperature analysis (Planck Collaboration XXIV 
2014) that error bars obtained with this simplified proce¬ 
dure are very close to optimal, provided the are pre¬ 
filtered with a simple diffusive inpainting technique (see 
Planck Collaboration XXIV 2014 for details). We find that this 
still holds true when we include polarization and pre-inpaint the 
T,Q,U input maps. Given its practical advantages in terms of 
speed and simplicity, we adopt this method in the following anal¬ 
ysis. 

A well-known, major issue with both Eq. (34) and (36) is 
that their direct implementation would require evaluation of all 
the bispectrum configurations from the data. The computational 
cost of this would scale like t)e totally prohibitive for 

high-resolution CMB experiments like Planck. The different bis¬ 
pectrum estimation techniques applied to our analysis are essen¬ 
tially defined by the approach adopted to circumvent this prob¬ 
lem. The advantage of having multiple independent implemen¬ 
tations of the optimal bispectrum estimator is twofold. Eirst, by 
cross-validating and comparing outputs of different pipelines, it 
strongly improves the robustness of the results. Second, different 
methods are complementary, in the sense that they have specific 
capabilities which go beyond simple /nl estimation. Eor exam¬ 
ple, the skew-Q method defined below facilitates the monitoring 
of NG foreground contamination, while the binned and modal 
estimators allow model-independent reconstruction of the data 
bispectrum, and so on. The skew-Q method enables the nature 
of any detected NG to be determined. Thus, the simultaneous 
application of all these techniques also allows us to increase the 
range and scope of our analysis. 
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In the following, we briefly outline the main features of 
the three optimal bispectrum estimation pipelines that are used 
for Planck measurements of /nl- We will only provide a 
short summary here, focused on the extension to polariza¬ 
tion data, referring the reader who is interested in more tech¬ 
nical aspects to our previous analysis of temperature data 
(Planck Collaboration XXIV 2014). 

3.1. KSW and skew-Ct estimators 

KSW and skew-C^ estimators (Komatsu et al. 2005; 
Munshi & Heavens 2010) can be used for bispectrum tem¬ 
plates that are written in factorizable (separable) form, i.e., 
as a linear combination of separate products of functions^ of 
ki, k 2 , and k^. This allows reduction of the three-dimensional 
integration over the bispectrum conflgurations into a product of 
three separate one-dimensional sums over ^i, ^ 2 , ^ 3 - This leads 
to a massive reduction in computational time (O(Vpix), where 
Vpix is the number of pixels in the map). The main difference 
between the KSW and skew-C^ pipelines is that the former 
estimates the /nl amplitude directly, whereas the latter initially 
estimates the so called “bispectrum-related power spectrum” 
(in short, “skew-C^”) function. Roughly speaking, the skew-C^ 
associates, with each angular wavenumber the contribution 
to the amplitude /nl (for each given shape) extracted from all 
triangles with one fixed side of size £. After resumming over the 
contributions from each ^-bin, the final point-like /nl estimate 
is obtained exactly as KSW. Equipping the KSW estimator with 
a skew-C^ extension can be particularly useful in the presence 
of (expected) spurious NG contaminants in the data. The slope 
of the skew-Q statistic is in fact shape-dependent and can be 
used to separate multiple NG components in the map. 


require separability of the starting theoretical shape in order to 
be applicable. Finally, after obtaining the data mode spectrum, it 
is possible to build a linear combination of the basis templates, 
using the measured amplitudes as coefficients, thus obtaining a 
model-independent full reconstruction of the bispectrum of the 
data. Of course the reconstructed bispectrum will be smoothed, 
as the estimator must use a finite number of basis templates. 

For this analysis, the modal method is implemented in 
two ways. One of them generalizes our previous temperature 
modal pipeline by expanding, for each shape, the correspond¬ 
ing TTT, EEE, TTE and EET bispectra. We then exploit sep¬ 
arability to build the covariance matrix of these expanded bis¬ 
pectra (Figuori et al. 2015), and to measure /nl efficiently using 
Eq. (36). This modal pipeline will be referred to throughout the 
paper as the “Modal 1” pipeline. 

The other implementation, which we will refer to as “Modal 
2”, utilizes a novel approach where the and are first or- 
thogonalized to produce new uncorrelated unit variance co¬ 
efficients. 



We then decompose the new bispectra as 


(41) 

(42) 


tXiX2X^ 



(43) 


3.2. Modal estimators 

Modal estimators (Fergusson et al. 2010a, 2012) are based on 
decomposing the bispectrum (both from theory and from data) 
into a sum of uncorrelated separable templates, forming a com¬ 
plete basis in bispectrum space, and measuring the amplitude 
of each. The evaluation of the amplitude for each template can 
be sped up by using a KSW approach (since the templates them¬ 
selves are separable by construction). All amplitudes form a vec¬ 
tor, also referred to as the “mode spectrum”. It is then possible 
to measure the correlation of the observed data mode spectrum 
with the theoretical mode spectra for different primordial shapes, 
in order to obtain estimates of the primordial /nl- Note also that 
the observed mode spectrum from data is theory-independent, 
and contains all the information from the data. Correlating the 
observed mode spectrum to theoretical mode vectors then al¬ 
lows the extraction of all the /nl amplitudes simultaneously. 
This makes modal estimators naturally suited for NG analyses, 
both when there are a large number of competing models to anal¬ 
yse, or when a model has free parameters through which we wish 
to scan (more than 500 shapes were analysed when applying this 
technique to Planck data). Another advantage is that by expand¬ 
ing into separable basis templates, the modal estimator does not 


^ We note that the local, equilateral, and orthogonal templates of 
Sect. 2 are separable. In fact, while the theoretical local NG models 
are manifestly separable, the equilateral and orthogonal templates of 
Eqs. (3) and (4) are factorizable approximations of the original non- 
separable shapes, that were derived exactly with the purpose of allow¬ 
ing the application of this type of estimator (Creminelli et al. 2006; 
Senatore et al. 2010). 


which can be constrained independently, since they are uncor¬ 
related. In this case the estimator then takes on a particularly 
simple form (Fergusson 2014). This new form is mathematically 
equivalent to the previous modal method, but involves signifi¬ 
cantly fewer terms in the estimator. However, due to the orthog- 
onalization procedure we cannot constrain the full EEE bispec¬ 
trum without further processing, just the additional part which is 
orthogonal to temperature. For this reason, although the “Modal 
2” T+E results incorporate all the polarization information, the 
EEE results alone are not presented here. 

In our analysis, both modal techniques (together with all the 
other estimators described in this section) were used to measure 
/nl for the three main shapes i.e., local, equilateral, and orthog¬ 
onal. Besides this, we optimized the two pipelines for different 
purposes. The “Modal 1” estimator was adopted to perform a 
large number of robustness tests of our results, especially in re¬ 
lation to the local, equilateral and orthogonal measurements. The 
“Modal 2” pipeline was instead mostly used to study a large 
number of “non-standard” primordial shapes (e.g., oscillatory 
bispectra). For this reason, each pipeline uses a different set of 
basis templates. The “Modal 1” estimator starts from a polyno¬ 
mial basis with 600 modes, and includes nine more modes that 
are the contributions from last scattering to the exact radial KSW 
expansion of the local, equilateral and orthogonal templates. The 
“Modal 2” expansion uses a high-resolution basis with 2000 
polynomial modes, augmented with a Sachs-Wolfe local bispec¬ 
trum template, in order to improve convergence efficiency in the 
squeezed limit. In this way, the high resolution estimator pro¬ 
vides the ability to scan across a wide variety of non-separable 
and oscillatory shapes, while the lower resolution pipeline gives 
efficient convergence in the /nl measurements for the standard 
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local, equilateral, and orthogonal shapes, offering rapid analy¬ 
sis for validation purposes. The “Modal 1” pipeline can also be 
generalized for the estimation of parity-odd bispectra, which is 
included in our analysis of non-standard shapes. 

3.3. Binned bispectrum estimator 

One can also use a binned estimator (Bucher et al. 2010, 2016), 
exploiting the fact that the theoretical bispectra of interest are 
generally smooth functions in ^-space. As a result, data and tem¬ 
plates can be binned in € with minimal loss of information, but 
with large computational gains from data compression. The data 
bispectrum in the binning grid is then computed and compared 
to the binned primordial shapes to obtain /nl- No KSW-like 
approach, which requires separability and mixing of theoreti¬ 
cal and observational bispectra in the computation, is required. 
Instead, the binned data bispectrum and the binned theoretical 
bispectrum and covariance are computed and stored completely 
independently, and only combined at the very last stage in a sum 
over the bins to obtain /nl- This means that it is very easy to 
test additional shapes or different cosmologies, and the data bis¬ 
pectrum can also be studied on its own in a non-parametric ap¬ 
proach. In particular the smoothed binned bispectrum approach, 
also used in this paper, investigates the (smoothed) binned bis¬ 
pectrum of the map divided by its expected standard deviation, to 
test if there is a significant bispectral NG of any type in the map. 
Another advantage of the binned bispectrum estimator is that the 
dependence of /nl on ^ can be investigated for free, simply by 
leaving out bins from the final sum. 

In more detail, the computation for the binned bispectrum 
estimator is based on Eqs. (38) and (39). However, instead of us¬ 
ing the reduced bispectrum , all expressions start from the 

alternative rotationally-invariant reduced bispectrum = 

' where h is defined in Eq. (40). The expression in 
Eq. (39) for the inner product remains the same when replac¬ 
ing h by B, except that the becomes h~^. The importance of 
is that it can be determined directly from maximally- 
filtered maps; 

(44) 

where 

M^in) = E (45) 

m 

Binning is then implemented by adding a sum over all i inside a 
bin to the expression for the filtered map given in Eq. (45), thus 
obtaining the observed binned bispectrum of the map 
with bin indices /i, h, h. The linear correction term is obtained 
in a similar way (and subtracted from the cubic term), but with 
two of the maps in Eq. (44) replaced by Gaussian simulations, 
and taking the average over a large number of those. 

The theoretical templates are binned simply by summing the 
exact expression over the ^s inside a bin,^ and the same is true for 

the covariance matrix The binning is op¬ 

timized in such a way as to maximize the overlap, defined using 
the inner product of Eq. (39) between the binned and the exact 
template for all shapes under consideration. Einally the estimate 

^ We note that the enormous computational gain of the binned bispec¬ 
trum estimator comes from the binned determination of the observed 
bispectrum; determining the theoretical bispectrum templates is fast, 
even when done exactly. 


of /nl is computed using Eq. (38), where the inner product now 
contains a sum over bin indices i instead of multipoles i, and 
the bispectra and covariance matrix are replaced by their binned 
versions. 

3.4. Data set and simulations 

In the following, we will apply our bispectrum estimation 
pipelines both to simulations and data, and consider a large 
number of shapes, either primordial or non-primordial in ori¬ 
gin. Simulations will be used for a wide range of purposes, from 
comparisons of the outcomes of different estimators, to tests of 
instrumental systematics and foreground contamination, as well 
as Monte Carlo evaluation of error bars. Eor this reason, many 
different sets of simulated maps will be used, with features that 
will vary, depending on the specific application, and will be de¬ 
scribed case by case throughout the paper. Most of the time, 
however, we will use the EEP8 simulation data set described 
in Planck Collaboration XII (2016), or mock data sets obtained 
by processing PPP8 maps in various ways. These are the most 
realistic simulations available, modelling the CMB sky and the 
instrumental effects of Planck to the best of our current knowl¬ 
edge. They have passed through the same steps of the component 
separation pipelines as the real sky map and are the same maps 
as used for the final validation of the estimators in Sect. 5.3. 

As far as actual data are concerned (Planck Collaboration I 
2016; Planck Collaboration II 2016; Planck Collaboration III 
2016; Planck Collaboration IV 2016; Planck Collaboration V 
2016; Planck Collaboration VI 2016; Planck Collaboration VII 
2016; Planck Collaboration VIII 2016) the maps analysed in 
this work are the Planck 2015 sky map, both in tem¬ 
perature and in E polarization, as cleaned with the four 
component separation methods SMICA, SEVEN, NILC, and 
Commander (Planck Collaboration IX 2016). As explained in 
Planck Collaboration VII (2016), the polarization map has had 
a high-pass filter applied to it, since the characterization of sys¬ 
tematics and foregrounds in low-^ polarization is not yet satis¬ 
factory. This filter removes the scales below t - lb completely, 
and those between I - lb and ^ = 40 partially. In all our analy¬ 
ses we use ^rnin = 40 for polarization, in order to be independent 
of the details of this filter. Eor temperature, we use ^min = 2. All 
the final cleaned maps are smoothed with a 5' Gaussian beam in 
temperature, and a 10' Gaussian beam in polarization. 

The maps are masked to remove the brightest parts of the 
Galaxy as well as significant point sources. The masks used 
are the common masks of the Planck 2015 release in tem¬ 
perature and polarization, which are the union of the con¬ 
fidence masks for the four component separation methods^ 
(Planck Collaboration IX 2016). The sky coverages are respec¬ 
tively /sky = 0.76 in temperature and /sky = 0.74 in polarization. 
The stability of our results as a function of the mask is investi¬ 
gated in Sect. 7.2, where we show that our temperature and joint 
temperature plus polarization results do not change significantly 
when we consider a larger sky coverage. 

In Sects. 7.1 and 7.3 we also compare the performance of 
different component separation methods, and conclude that, with 


^ We note that the Planck collaboration produced two slightly dif¬ 
ferent sets of union masks (see Planck Collaboration IX 2016 for de¬ 
tails). We choose to adopt the more conservative set in this paper, as 
we found that the agreement between different component separation 
methods significantly increases with these masks when we measure /^l 
of shapes that peak in the squeezed limit (while the differences are very 
small in other cases). 
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respect to bispectrum estimation, the most accurate results are 
obtained using the SMICA data set. As already done for the 2013 
release, we will thus consider SMICA as our main data set, using 
the other methods for important cross-checking purposes. 

If we consider only temperature, current SMICA data become 
noise dominated at ^ ^ 2000, while previous nominal mis¬ 
sion data were noise dominated at ^ ^ 1700. The mask used 
for the 2013 release was also slightly larger than the current 
one (/sky = 0.73 in 2013 vs. /sky = 0.76 in 2015). Since the 
/nl signal-to-noise ratio, as quantified by the Fisher Information 
Matrix, scales as (S/N) cx i in the signal dominated regime 
and saturates in the noise-dominated regime, we expect an im¬ 
provement in our /nl temperature constraints of about 20% 
when going from the 2013 nominal mission release to the current 
results. Adding polarization and accounting for all possible TTE, 
EET, and EEE bispectra produces further improvements. Since 
we are neglecting the first 40 polarization multipoles, such im¬ 
provements are expected to be fairly small for shapes peaking 
in the squeezed limit, and more pronounced for equilateral type 
bispectra. A Fisher matrix approach shows that error bars are ex¬ 
pected to improve by about 10 % for the local shape and about 
40 % for the equilateral shape. This is in good agreement with 
our actual measurements, as can be seen from the results pre¬ 
sented in Sect. 6 onwards. 

3.5. Data analysis settings 

Now we detail the general setup adopted for the analysis of 
Planck 2015 data by the four different optimal bispectrum es¬ 
timation pipelines, described in previous sections. 

As already explained, inpainting of the masked regions of 
the sky is a preliminary data filtering operation that all pipelines 
must perform, in order to retain optimality. We found that the 
inpainting method used in 2013 (Planck Collaboration XXIV 
2014) for temperature maps still works well when polarization 
is included (note that it is the original T, 2, and U maps that 
should be inpainted, not the derived E map). We adopt a simple 
diffusive inpainting method. First the masked regions of the map 
are filled with the average value of the rest of the map. Then 
the value of each masked pixel is replaced by the average value 
of its (generally eight) direct neighbour pixels. The latter step is 
repeated a fixed number of times (2000).^^ Relevant, final com¬ 
putations in map space (see e.g., Eq. 44) are always done after 
remasking, so that the inpainted regions of the map are not used 
directly. The relevance of the inpainting procedure is that it re¬ 
duces the effect of the sharp edges and the lack of large-scale 
power inside the mask leaking into the rest of the map during 
harmonic transforms. 

For the linear correction term and to determine error bars, we 
use the FFP8 simulations (see Planck Collaboration XII 2016, 
and Sect. 3.4), filtered through the different component separa¬ 
tion pipelines, using the same weights as used for the actual data 
when co-adding frequency channels. To compute all theoretical 
quantities (like the bispectrum templates and the ISW-lensing 
bias) we use the Planck 2015 best-fit cosmological parameters 


For bispectrum purposes we found no difference between the results 
when performing the procedure without a buffer (the so-called “Gauss- 
Seidel” method, where amongst the neighbours will be pixels both at 
the current and at the previous iteration) and with a buffer (the so-called 
“Jacobi” method, where all neighbour pixels will be at the previous iter¬ 
ation), except that the former converges faster. We found 2000 iterations 
to work well in the ‘Gauss-Seidel’ case. 


as our fiducial cosmology. However, results are quite insensitive 
to small changes in these parameters. 

As pointed out in Sect. 3.4, low-^ multipoles are filtered out 
of the input polarization data set, so that all of our analyses will 
use 7’min = 40 in polarization, and ^min = 2 in temperature. The 
choice of ^max is dictated by the angular resolution of the cleaned 
maps, which is 5' in temperature, and 10' in polarization, and 
by the fact that the temperature data become noise-dominated at 
i ^ 2000, while the polarization information saturates around 
I ^ 1000. The KSW and binned estimators use ^max = 2500 
for temperature, while the modal estimators use ^max = 2000. 
As shown explicitly in Sect. 7.4, results are completely stable 
between I = 2000 and i = 2500, so that this has no impact on 
/nl- Similarly the binned estimator uses T’max = 2000 for polar¬ 
ization, while the other estimators use ^max = 1500, but again 
Sect. 7.4 shows that this difference is unimportant. The estima¬ 
tors also differ in the number of maps used to compute the linear 
correction term and the error bars, but generally it is of the or¬ 
der of 200. This difference is due to the different convergence 
properties of the estimators, some converging faster than others. 

The binned bispectrum estimator uses 57 bins^^ for the anal¬ 
ysis, which were determined by optimizing the correlation be¬ 
tween the exact and the binned templates for the different shapes 
in temperature and polarization, as well as the full combined 
case. This is equivalent to minimizing the variance of the dif¬ 
ferent /nl parameters, where we focus in particular on the pri¬ 
mordial shapes. 

As previously explained, we use two different versions of 
the polarized modal pipelines, called “Modal 1” and “Modal 2” 
in the paper. Besides technical and conceptual implementation 
differences, the two modal estimators also use different sets of 
basis templates. The “Modal 1” pipeline uses 600 polynomial 
modes, plus nine “KSW radial modes”, computed at last scat¬ 
tering, while “Modal 2” has a basis formed by 2000 polynomial 
modes, augmented with a Sachs-Wolfe local bispectrum tem¬ 
plate. Due to the way polarization is implemented in the “Modal 
2” pipeline, it cannot determine results for £'-only. More de¬ 
tails and explanations of the different choices are provided in 
Sect. 3.2. 

As already stressed, the use of several independent bispec¬ 
trum estimators, and several completely independent component 
separation methods allows a remarkable level of cross-validation 
of our results in order to establish their robustness. The fact that 
the bispectrum estimators are statistically equivalent and pro¬ 
duce practically optimal results will be established in Sect. 5. 
The validation of the component separation methods is described 
in Planck Collaboration IX (2016) and Sect 7. 


4. Non-primordial contributions to the CMB 
bispectrum 

Here we investigate several bispectra of non-primordial origin 
that are expected to be present in the data, and quantify their im¬ 
pact on our /nl results. We devote particular attention to assess¬ 
ing potential biases that these NG signals might induce on the 


The boundary values of the bins are: 2, 4, 10,18, 30, 40, 53, 71, 99, 
126, 154, 211, 243, 281, 309, 343, 378, 420, 445, 476, 518, 549, 591, 
619, 659, 700, 742, 771, 800, 849, 899, 931, 966, 1001, 1035, 1092, 
1150, 1184, 1230, 1257, 1291, 1346, 1400, 1460, 1501, 1520, 1540, 
1575,1610,1665,1725,1795,1846,1897, 2001, 2091, 2240, and 2500 
(i.e., the first bin is [2, 3], the second [4, 9], etc., while the last one is 
[2240, 2500]). 
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primordial bispectra. When forecasting such biases, we assume Table 1. Bias in the three primordial /nl parameters due to the 
the data analysis settings discussed in Sects. 3.4 and 3.5. lensing-ISW signal for the four component separation methods. 


4.1. Non-Gaussianity from the lensing-ISW bispectrum 

The correlation between the gravitational lensing of the CMB 
anisotropies and the integrated Sachs-Wolfe (ISW) effect gives 
rise to a secondary CMB bispectrum — characterized by 
an oscillatory behaviour and peaked on squeezed configura¬ 
tions — that is a well-known contaminant to the primordial 
NG signal (Hanson & Lewis 2009; Mangilli & Verde 2009; 
Lewis etal. 2011; Mangilli et al. 2013). The temperature- 
only 2013 Planck results (Planck Collaboration XXIV 2014; 
Planck Collaboration XIX 2014; Planck Collaboration XVII 
2014) showed evidence for the first time for the lensing- 
ISW CMB bispectrum and associated bias. Based on the 
same methodology used for the 2013 Planck data analysis 
(Planck Collaboration XXIV 2014), here we update the compu¬ 
tation of the lensing-ISW bispectrum and its bias to include the 
full mission temperature and polarization data. 

As shown by Cooray & Melchiorri (2006), the direct 
lensing-ISW correlation in ^’-polarization due to rescattering of 
the temperature quadrupole generated by the ISW effect is neg¬ 
ligible. However, as explained in Lewis et al. (2011), there is 
an important correlation between the lensing potential and the 
large-scale £’-polarization generated by scattering at reioniza¬ 
tion. Because the lensing potential is highly correlated with the 
ISW signal, this also leads to a non-zero lensing-ISW bispec¬ 
trum in polarization. 

To determine the amplitude parameter of the lensing- 

ISW bispectrum, one simply inserts the theoretical template for 
this shape into the general expression of Eq. (38). The template 
is given by (Hu 2000; Lewis et al. 2011) 


T XiX2X3,LISW_ rX\ ^Xj,^^X\X 2 ^X\ 

rX2 -X 30 ^XiX 2 rX2 

^Xi0^X2X3 rX^ ^X20^XiX3 rX^ 


(46) 


Here and are the temperature/polarization-lensing po¬ 
tential cross power spectra, and the tilde on and 

indicates that it is the lensed TT, TE, or EE power spectrum. 
The functions are defined by 


fE 


- [^2(^2 + 1 ) + ^3(^3 + 1 ) “ A(A + 1 )] , 

- [^2(^2 + 1 ) + ^3(^3 + 1 ) “ A(A + 1 )] 


/G 4 

\2 b -2 


ix h ^3 V 

0 0 0 / 


(47) 


if -I- ^2 + ^3 is even and , ^ 2 , ^3 satisfy the triangle inequality, 
and zero otherwise. 

In this paper our main concern with the lensing-ISW bispec¬ 
trum is not so much to determine its amplitude (although that is 
also of great interest), but to compute its influence on the primor¬ 
dial shapes. The bias due to the lensing-ISW bispectrum on 
the estimation of a given primordial amplitude given by 


where the inner product is defined in Eq. (39). 


(48) 


lensing-ISW /nl bias 


Shape SMICA SEVEM NILC Commander 


r Local . 7.5 7.5 7.3 7.0 

r Equilateral. 1.1 1.2 1.3 1.8 

P Orthogonal .... -27 -27 -26 -26 


E Local . 1.0 1.1 1.0 1.1 

E Equilateral. 2.6 2.7 2.5 2.9 

E Orthogonal .... -1.3 -1.3 -1.2 -1.5 


r+E Local . 5.2 5.5 5.1 4.9 

E+E Equilateral .. 3.4 3.4 3.4 3.6 

E+E Orthogonal .. -10 -11 -10 -10 


Table 2. Results for the amplitude of the lensing-ISW bis¬ 
pectrum from the SMICA, SEVEM, NILC, and Commander 
foreground-cleaned maps, for different bispectrum estimators. 
Error bars are 68 % CL; see the main text for how they have 
been determined. 




lensing-ISW amplitude 


Method 

SMICA 

SEVEM 

NILC 

Commander 

E 

KSW . . . . 
Binned .. . 
Modal2 . . 

0.79 + 0.28 
0.59 + 0.33 
0.72 + 0.26 

0.78 + 0.28 
0.60 + 0.33 
0.73 + 0.26 

0.78 + 0.28 
0.68 + 0.33 
0.73 + 0.26 

0.84 + 0.28 
0.65 + 0.36 
0.78 + 0.27 

E+E 

Binned .. . 

0.82 + 0.27 

0.75 + 0.28 

0.85 + 0.26 

0.84 + 0.27 


The values for the bias are given in Table 1. It should 
be noted that these are the results as computed exactly with 
Eq. (48). They can differ slightly from the ones used in e.g.. 
Table 10, where each estimator adopts values computed us¬ 
ing the approximations appropriate to the method. However, 
these differences are completely insignificant. As seen already in 
Planck Collaboration XXIV (2014), for E-only the bias is very 
significant for local and to a lesser extent for orthogonal NG. Eor 
local NG the bias is larger than the error bars on /nl- We see that 
for E-only the effect is non-zero but not significant. Eor the full 
E-rE case, the bias is smaller than for E-only, but large enough 
that it is important to take into account. 

The results for can be found in Table 2. The polarized 
version of the template has only been implemented in the binned 
bispectrum estimator. Error bars have been determined based on 
EEP8 simulations as usual. The KSW estimator implements 
the lensing-ISW template exactly, while the binned and modal 
estimators use approximations, as explained in Sect. 3. In partic¬ 
ular for the binned estimator the correlation between the binned 


The average value of the lensing-ISW amplitude determined from 
the FFP8 simulations is around 0.85 of the expected value. This value 
is very consistent across bispectrum estimators and component separa¬ 
tion methods, which provides a useful consistency test in its own right. 
Except for this effect, all other tests on the temperature FFP8 maps show 
them to be very robust and to behave as expected, for example in the 
determination of the lensing-ISW bias on the local shape. We took this 
effect into account by increasing all error bars in the table by the appro¬ 
priate factor (i.e., dividing them by ^ 0.85). 
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Fig. 1. The skew-Q spectrum for the lensing-ISW effect (red 
line with data points), from the temperature map. The blue curve 
is the theoretically-expected spectrum. Note that the points be¬ 
yond I = 1500 are significantly correlated. 



Fig. 2. The skew-Q spectrum for unclustered point sources (red 
line with data points), from the temperature map. The blue curve 
is the theoretical spectrum, given the amplitude determined with 
the KSW estimator. 


and exact lensing-ISW template is relatively low, since it is a dif¬ 
ficult template to bin (unlike all the other templates considered 
in this paper), which is reflected in the larger error bars. Tests 
performed on FFP8, as well as other tests, demonstrate that the 
lower correlation does not lead to a bias compared to the other 
estimators. We will use the KSW results to draw our conclu¬ 
sions. 

We see that temperature results from the full mis¬ 
sion are consistent with the 2013 nominal mission 
(Planck Collaboration XXIV 2014). Including polarization 
yields results that also appear consistent and decrease the 
error bars. However, for now the T+E conclusions should be 
considered preliminary, for the reasons related to polarization 
data discussed in detail in Sects. 6 and 7. The error bars will also 
improve when measured with the other bispectrum estimators. 
As already seen in 2013, the values for are slightly low 

compared to the expected value of 1, but not significantly so. 
On the other hand, the detection of the lensing-ISW bispectrum 
is significant, even with our conservative rescaling of the error 
bars. The hypothesis of having no lensing-ISW bispectrum 
is excluded at 2.8 cr using temperature alone, and improves 
to 3.0 cr with the current preliminary result when including 
polarization. As mentioned above, the latter result is likely to 
improve with further analysis of the Planck data. In Fig. 1 we 
present the results of the skew-C^ analysis for lensing-ISW NG 
for the T map, which illustrates that the instrument and data 
processing are not removing this expected NG signal from the 
data. 

4.2. Non-Gaussianity from extragalactic point sources 

The auto-bispectra of extragalactic point sources are a potential 
contaminant to primordial NG estimates at Planck frequencies. 
The basic modelling and methodology of this section follows the 
corresponding section in Planck Collaboration XXIV (2014). 

Extragalactic point sources are divided into populations 
of unclustered and clustered sources. The former are radio 
and late-type infrared galaxies (see e.g., Toffolatti et al. 1998; 
Gonzalez-Nuevo et al. 2005), while the latter are dusty star¬ 


forming galaxies constituting the cosmic infrared background 
(GIB; Lagache et al. 2005). The contamination due to both types 
of sources in NG estimators is handled via dedicated bispectrum 
templates which are fitted jointly with the primordial NG tem¬ 
plates. 

The unclustered sources have a white noise distribution, and 
hence constant polyspectra. Their reduced angular bispectrum 
template is thus 

= constant. (49) 

This constant is usually noted Z^ps or Z?src m the literature (e.g., 
Komatsu & Spergel 2001). This constant template is valid in po¬ 
larization as well as temperature, since the polarization angles of 
point sources are less clustered than the source density. However, 
since not all these point sources are polarized, we do not measure 
the same sources in temperature and in polarization. In fact, there 
is no detection of the bispectrum of unclustered point sources in 
the cleaned Planck polarization map, unlike in the temperature 
map, where Table 3 (binned bispectrum estimator) and Fig. 2 
(skew-Qs) show a clear detection. 

The clustered sources (GIB) have a more complex bispec¬ 
trum in temperature, refiecting the distribution of the large- 
scale structure and the clustering of galaxies in dark matter 
halos (Argueso et al. 2003; Lacasaetal. 2012; Crawford et al. 
2014). The Planck results have allowed the measurement of the 
GIB bispectrum at frequencies 217, 353, 545 GHz in the range 
i ^ 200 - 700 (Planck Collaboration XXX 2014). In this multi¬ 
pole range, a power law was found to fit the measurement, with 
an exponent consistent between frequencies. However, at lower 
multipoles theoretical models for the GIB power spectrum (e.g., 
Planck Collaboration XXX 2014) and bispectrum (Lacasa et al. 
2014; Penin et al. 2014) predict a fiattening of the GIB power. 
We thus take the TTT GIB bispectrum template to be a broken 
power law. 


b 


CIB 


oc 


(1-1- 7'i/4reak)(l + ^2/^break)(l + ^s/^break) ^ 

(1 + ^o/^break)^ 


(50) 


where the index is ^ = 0.85, the break is located at 4reak = 
70, and = 320 is the pivot scale for normalization. Dusty 
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Table 3. Joint estimates of the bispectrum amplitudes of unclus¬ 
tered and clustered point sources in the cleaned Planck tempera¬ 
ture map, determined with the binned bispectrum estimator. The 
error bars have been determined using FFP8 simulations. 


Map 

^s/dO-"'') 

Acib/(10 ^^) 

SMICA . 

5.6 + 2.7 

0.4 + 1.4 

SEVEM . 

7.9 + 2.8 

0.8 + 1.4 

NILC. 

9.3 + 2.7 

-0.3 + 1.4 

Commander . . . . 

5.9 + 3.2 

1.4 + 1.6 


star-forming galaxies emit with a low polarization fraction, and 
polarization correlates only over the smallest scales, so that the 
CIB is negligibly polarized. We thus take vanishing templates 
for its polarized bispectra 


iCmjTE _ iCmjEE _ iCm,EEE _ ^ 

^ p p p ^ p p p ^ p p p 

€1C2^3 ^1^2^3 


(51) 


Both point source templates, Eqs. (49) and (50), have been 
implemented in the binned bispectrum estimator described in 
Sect. 3. The results for these two templates applied to the Planck 
temperature map, cleaned with the four component separation 
methods, can be found in Table 3. Since the two templates are 
highly correlated, the results have been determined in a com¬ 
bined analysis. The results have also been determined jointly 
with the primordial local, equilateral, and orthogonal templates, 
and the lensing-ISW bias has been subtracted, but all of this 
makes a negligible difference. Contamination from unclustered 
sources is detected in all component-separated maps. However, 
Acib is not detected. 

The order of magnitude of the bispectrum amplitudes found 
in Table 3 is consistent with expectations. Indeed, for radio 
sources at 217 GHz and with a flux cut based on the Planck 
ERCSC (Planck Collaboration VII 2011), Lacasa & Aghanim 
(2014), the forecasted Z^ps is approximately 2 x 10“^^. For the 
CIB, the Planck 2013 measurement (Planck Collaboration XXX 
2014) at 217 GHz gives Acm - 6x10“^^ when translated into di¬ 
mensionless units. The results reported in Table 3 are consistent 
at the order-of-magnitude level with these estimates, although 
they are lower, because we are analysing cleaned maps. 

The unclustered point source and CIB templates are highly 
correlated, at 93 %. For this reason it was not deemed a priority 
for the other bispectrum estimators to implement the CIB tem¬ 
plate as well. Moreover, both point source templates are negligi¬ 
bly correlated with the primordial NG templates and the lensing- 
ISW template (the maximum being the correlation between equi¬ 
lateral and CIB templates at 2.7 %, while correlations with the 
unclustered point source template are well below 1 %). For this 
reason, and despite the detection of point sources in the cleaned 
maps, it makes no difference for the primordial results if point 
sources are included in a joint analysis or completely neglected. 

An additional contaminant to the cosmological CMB bis¬ 
pectrum arises from the correlation between the gravitational 
lensing of the CMB anisotropies and the CIB anisotropies. This 
correlation was detected in Planck Collaboration XVIII (2014) 
using an optimal cross-spectrum estimator. The CIB-lensing 
bispectrum might couple with any of the primordial shapes. 
However, the amplitude of the CIB bispectrum is predicted to 
be small in the cleaned Planck maps and it has actually not been 
detected (see Table 3). The CIB-lensing bispectrum signal is fre¬ 
quency dependent, and it is mostly dominant in the very high 
Planck frequencies (see e.g., Curto et al. 2014). 


4.3. Non-Gaussianity from residuals of the deglitching 
processing 

Cosmic rays interacting with the cryogenic detectors induce 
spikes in timelines. These high-amplitude, fast-rising signals are 
followed by a decay tail. We observe three families of glitches, 
characterized by their temporal shape. The amplitude and time 
constants of the decays depend on which part of the satellite is hit 
(Catalano et al. 2014; Planck Collaboration X 2014). These ran¬ 
dom events are Poisson-distributed in time and produce highly 
non-Gaussian systematics. 

A method has been developed to remove them directly at the 
time-ordered information (TOI) level. This process is performed 
iteratively, and is described in detail in Planck Collaboration X 
(2014). The short glitches are just flagged in the data, whereas 
for the long ones only the fast part is flagged, and the long tail is 
substracted from the timeline. This procedure is not perfect, and 
there are residuals from the potentially biased errors in the fit, 
and the undetected glitches under the threshold of 3.2(T of the 
TOI noise rms. They could in principle produce a non-Gaussian 
signal in the final map. In addition, these residuals could interact 
with the mapmaking procedure at the destriping level, since the 
error on the offset determination could be non-Gaussian due to 
undetected glitches or a possible bias in the errors of the removal 
of tails. This is important, because in more than 95 % of the TOI 
data, tails have been subtracted. 

To estimate the effect of these residuals on the determi¬ 
nation of NG, we created two sets of simulations (one in¬ 
cluding glitches and the other not) for every bolometer of the 
143 GHz channel. We generated Gaussian CMB maps, and 
applied the full TOI processing with a realistic instrumental 
noise (Planck Collaboration VII 2016). In the simulations with 
glitches, we added glitches at the TOI level, following the prop¬ 
erties measured in the data, and cleaned them with the procedure 
applied to the data. For the simulations without glitches, we have 
the same CMB and noise realization, but no glitches added at the 
TOI level. 

We estimated the bias caused by glitches on the measurement 
of /nl using the binned bispectrum estimator. The bias on /nl 
induced by the glitch residuals g on a map T, including noise and 
CMB is given by + g) - fm.(T)), where the noise in the 

weighting of the estimator is determined from the simulations 
with glitches (as it would be for the data). Results are shown 
in Table 4. For most shapes, we detect no significant bias. The 
higher signal and high dispersion for the local shape might be 
due to a mis-calibration of the linear correction. In any case, for 
all shapes the bias due to glitches is a negligible correction to 
the value of /nl, given its error bars, and we will not take it into 
account in the remainder of the paper. 


5. Validation tests 

During the work for the 2013 release, culminating in the NG 
results of Planck Collaboration XXIV (2014), the advantage of 
having multiple independent bispectrum estimator implemen¬ 
tations was amply demonstrated. This allows for very useful 
cross-checking of results, both during development and for the 
final analysis, thus greatly improving the robustness of and con¬ 
fidence in the final results. For this new release we followed 
the same procedure, with the same three principal bispectrum 
estimators: KSW; binned; and modal, all of which had their 
pipelines updated to handle polarization data in addition to tem¬ 
perature. 
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Table 4. Results on the impact of cosmic ray residuals on the estimation of /nl at 143 GHz, determined using the binned bispectrum 
estimator. We produced 10 simulations. We report the mean of the bias defined in the text, and the error on this mean. We also show 
the Fisher error bars on /nl for these simulations. 


Local Equilateral Orthogonal Diffuse PS (x 10^^) Lens-ISW 


T -only 

bias mean. 1.1 + 0.6 0.8 + 1.6 -1.0+ 0.7 0.5 +0.2 0.01 + 0.01 

o-f^^ . 5.2 64 34 2 0.2 


F-only 

bias mean. 2.4 + 5.8 -9.6 + 8.7 -7.1 + 14.8 0.0 +0.1 -3.0 +1.4 

. 38 157 90 0.6 7.8 


T+E 

bias mean. 1.8 + 1.0 -5.1 + 2.2 0.1 + 1.5 0.01 + 0.04 0.01 + 0.01 

o-f^^ . 4.4 43 22 0.3 0.2 


Beyond the usefulness of cross-checking, the three esti¬ 
mators complement each other and have different strengths. 
The KSW estimator can treat separable bispectrum templates 
without approximation, but it is more work to add new tem¬ 
plates and non-separable templates cannot be handled at all. The 
binned and modal estimators can reconstruct the full bispectrum 
(smoothed in different domains), while the skew-Q extension 
of the KSW estimator can be used to investigate the bispectrum 
beyond /nl- The binned bispectrum estimator is the fastest on a 
single map or a set of unrelated maps, but becomes slower than 
the other two on a large set of realizations based on the same set¬ 
tings, because the linear correction term cannot be precomputed. 
The modal estimator can investigate a wide selection of oscil¬ 
lating or otherwise rapidly changing bispectrum templates that 
would be difficult to bin, while the binned bispectrum estimator 
can quickly implement and determine the /nl of an additional 
template or the effect of a different cosmology if the binned bis¬ 
pectrum of the maps has already been computed. The binned es¬ 
timator gets the dependence of /nl on I for free with its results, 
while the modal estimator allows for a statistical investigation of 
the mode coefficients. 

In this Section we show some of the validation tests, in par¬ 
ticular for polarization. In Sect. 5.1 we investigate the agreement 
between estimators, map-by-map, on sets of successively more 
realistic maps. In Sect. 5.2 we show that the estimators are unbi¬ 
ased in the presence of a non-zero /nl- Finally, in Sect. 5.3 we 
show that the estimators are essentially optimal on a set of the 
most realistic Planck simulations available, which are those used 
to compute the error bars on our final results. 

5.1. Agreement between estimators on a map-by-map basis 

The maps used in this subsection are realistic simulations of the 
CMB (at resolution A^side = 2048) but without any foregrounds. 
They do not contain any primordial NG, but do include ISW- 
lensing. Since the final FFP8 simulations were not yet available, 
the main goal was to make sure that the estimators agreed with 
each other, not only on average, but also on a map-by-map ba¬ 
sis. For this purpose it was enough to look at only 49 maps. 
Establishing optimality of the estimators requires a larger num¬ 
ber of maps, and is shown on the FFP8 simulations in Sect. 5.3. 

In our first test we include the effect of the 143 GHz beam, 
but in other respects the simulations are ideal (no noise, and no 
mask). The analysis used T’max = 2000 for both T and E. The 
results for the average over the maps for the KSW, binned, and 


both modal estimators, as well as for the difference between each 
estimator and KSW, are shown in Table 5. The shapes are as¬ 
sumed to be independent in this analysis, which means that the 
bias on the local shape due to the ISW-lensing effect is clearly 
visible. Results are shown for T-only, £’-only, and the full com¬ 
bined T-\-E analysis. Note that the second modal implementation 
cannot compute results for E alone. One clearly sees that the re¬ 
sults agree very well. It is also interesting to note that in this 
ideal noiseless case, one can actually determine /nl more accu¬ 
rately from polarization alone than from temperature alone. This 
is due to the narrower transfer function in polarization, so that 
the primordial bispectrum is less smoothed in its projection to 
two-dimensional harmonic space. 

The second test is identical to the first, except that we add re¬ 
alistic anisotropic noise realizations to the full-sky maps, based 
on the 143 GHz channel. The estimators now require the use of 
the linear correction term, and results are shown in Table 6. The 
agreement is still very good, although slightly worse than in the 
ideal case, as expected. The fact that the error bars for the T -only 
local case here are actually a bit smaller than in the ideal case is 
an artefact of the small number of maps; i.e., the error bars have 
not completely converged yet. On the other hand, the fact that 
the error bars for £’-only are much larger than in the ideal case is 
a real effect; the Planck single-frequency polarization maps are 
noise-dominated. 

Finally, the third test is identical to the second, except that 
we now also add a mask. The mask chosen is realistic, based on 
the union of the confidence masks provided by the SMICA, NILC, 
SEVEN, and Commander methods for this particular set of simu¬ 
lations. It contains both a Galactic and a point source part. The 
temperature mask leaves 79 % of the sky unmasked, while the 
polarization mask leaves 76 %. The results are shown in Table 7, 
while the map-by-map comparison is given in Fig. 3. From the 
table we see that the agreement between the different bispectrum 
estimators is still very good and only slightly degraded when 
compared to the previous case. The typical discrepancy between 
the bispectrum estimators, even in this most realistic case, is less 
than about a third of the uncertainty on /nl- This is apparent in 
the map-by-map comparison of Fig. 3. 

5.2. Validation of estimators in the presence of primordial 
non-Gaussianity 

After the map-by-map comparison of the previous section, we 
next want to make sure that the estimators are unbiased. For 
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Table 5. Results from the different /nl estimators for the set of CMB simulations described in Sect. 5.1 in the ideal case without 
noise or mask. Both the results for the estimators individually and for the differences with KSW are given, for T-only, £’-only, and 
the full combined T+E analysis. The shapes are assumed independent. 


/nl 


Shape 

KSW 

Binned 

Modal 1 

Modal 2 

B - KSW Ml - KSW 

M2 - KSW 

T Local . 

T Equilateral. 

T Orthogonal . 

_ 7.6 + 5.4 

.... 7 +53 

_ -22+27 

7.4+ 5.6 

5 +58 
-22 + 28 

7.4+ 5.1 

6 +53 
-22 + 27 

7.2+ 5.7 

8 +56 

-17 +30 

-0.3+ 0.6 
-2 + 12 
0.5+ 9.4 

-0.2 + 0.4 
-1.0+ 8.4 
-0.2 + 4.2 

-0.5+ 2.2 
0+17 
5+11 

E Local . 

F Equilateral. 

E Orthogonal . 

_ -0.9 + 4.1 

_ -9 +42 

. . . . 4 +13 

-1.3+ 3.4 
-10 +42 
5+13 

-0.9+ 3.7 
-10 +40 
4+12 


-0.3+ 2.9 
-1 + 11 

0.1 + 3.8 

0.1 +0.5 
-0.7 + 9.2 
-0.3 + 2.7 


T+E Local . 

T+E Equilateral .... 
T+E Orthogonal .... 

_ 2.2+ 3.1 

.... 0 +20 
_ -4 +10 

1.5+ 2.5 

2 +22 
-4+9 

2.1+ 2.8 

3 +21 

-6+9 

2.0+ 3.3 

0 +23 
-5 + 12 

-0.6+ 1.0 
1.4+ 5.8 
0.3+ 2.2 

0.0 + 0.8 

2.3 + 7.3 
-1.1 + 3.1 

-0.2+ 1.9 
0+12 
-1.0+ 7.1 

Table 6. As Table 5, but with noise and no mask. 





/nl 




Shape 

KSW 

Binned 

Modal 1 

Modal 2 

B - KSW 

Ml - KSW 

M2 - KSW 

T Local . 

r Equilateral. 

T Orthogonal . 

6.7 + 4.8 

11+61 
. . -19 + 31 

6.4 + 5.2 

12 + 65 
-18 + 34 

6.7 + 4.7 

9+63 
-20 + 32 

7.0+ 5.3 

12 + 62 
-18 +35 

-0.3+ 1.0 
1 + 15 

1 + 12 

0.1+ 0.4 
-1.9+ 9.6 
-1.3+ 5.1 

0.3+ 1.2 

1 + 12 
0.8+ 8.8 

E Local . 

F Equilateral. 

E Orthogonal . 

-2+29 
1+191 
-6 + 101 

-4 + 28 
-18 + 195 

0 + 107 

-1 + 29 

-6 + 200 
-6 +102 


-2 + 12 
-19 +47 

6 +25 

0.4+ 5.5 
-7 +23 
-0.3 + 10 


T+E Local . 

T+E Equilateral . 

T+E Orthogonal. 

4.9 + 4.2 

13+46 
. . -11 + 22 

4.5 + 4.4 

11 + 49 

-11 +24 

5.0 + 4.2 

9+48 
-13 + 22 

4.9+ 4.9 

13 +47 

-11 +24 

-0.4+ 1.2 
-2 + 10 
0.0+ 7.3 

0.1+ 1.5 
-4 + 13 

-1.3+ 7.1 

-0.0+ 1.2 
-0.3 + 7.0 
0.7 + 4.5 

Table 7. As Table 5, but with noise and a mask. 





/nl 




Shape 

KSW 

Binned 

Modal 1 

Modal 2 

B - KSW 

Ml - KSW 

M2 - KSW 

T Local . 

T Equilateral. 

T Orthogonal . 

6.5+ 5.1 

11+73 
. . -22 + 37 

6.1 + 5.3 

9+75 
-21 + 37 

6.4 + 5.0 

6+76 
-23 + 36 

6.0+ 5.3 

11 +70 

-20 + 37 

-0.4+ 1.5 
-2 + 19 
2+14 

-0.1+ 0.7 
-5 + 14 

-0.9+ 6.1 

-0.5+ 1.3 
0+12 
2.6+ 9.2 

E Local . 

F Equilateral. 

E Orthogonal . 

4+36 
. . -32 + 242 

-9 + 138 

0+35 
-49 + 209 
-7 + 139 

5+37 
-38 + 246 
-7 + 142 


-4 + 16 

-17 +88 

2 +45 

1 + 13 

-6 +34 
2+19 


T +E Local . 

T +E Equilateral . 

T +E Orthogonal. 

5.1 + 5.3 

19+50 
. . -12 + 25 

4.2+ 5.1 

16 + 50 
-11 + 26 

4.8 + 5.0 

15 + 53 
-13 + 25 

4.5+ 5.2 

16 +45 
-11 +23 

-1.0+ 1.7 
-3 + 14 

1.9+ 8.7 

-0.3+ 1.7 
-4 + 19 
-1.0+ 9.9 

-0.6+ 1.3 
-3.2+ 9.8 
1.4+ 5.9 


this purpose we prepared a different set of 100 T and E CMB 
simulations, still with cosmological parameters as determined 
by Planck. This time ISW-lensing is not present, but there is 
a nonzero local /nl = 12. To these maps we add the same 
beam, anisotropic noise, and mask as before. We again take 
^max = 2000, and the results are given in Table 8. 

We see that all the estimators correctly recover the input 
value, both in temperature and in polarization. The results for 
the equilateral and orthogonal shapes are consistent with the fact 
that those templates have a non-zero correlation with the local 


shape (the Table gives the results for an analysis where all shapes 
are assumed independent). For example, a joint analysis of the 
T-rE’binned estimator gives= 11.5+6.4,/^^^^^ = -7.5+51, 
and /nl^° = -0.4 ± 29. Except for the first modal estimator in 
£’-only (due to an insufficient number of maps in the linear cor¬ 
rection term), we also find that the error bars for the bispectrum- 
based estimators are very close to the Fisher errors. Note that a 
slight increase in the error bars compared to Fisher estimates is 
expected for the local shape in T -only and in T +£", due to the 
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Fig. 3. Map-by-map comparison of the results from the different estimators for local (left), equilateral (centre), and orthogonal 
(right) /nl (taking the shapes to be independent), for the third set of simulations described in Sect. 5.1, including both noise and a 
mask. Results are shown for T-only (top), £’-only (centre), and the full combined T+E case (bottom). The legend for the estimators 
can be found in the top right figure. The horizontal solid line is the average value of all maps for KSW, and the dashed and dotted 
horizontal lines correspond to ±lcr and ±2cr deviations, respectively. 


signal being significantly different from zero there (the Fisher 
error bars for the local case are 5.8 for T-only, 26 for £’-only, 
and 5.0 for T-\-E). Hence the estimators are effectively optimal, 
as will be illustrated in more detail in the next section. 


5.3. Validation of estimators on realistic Planck simulations 

As a final validation test, we ran our estimators on a large set of 
the most realistic simulations available. These are the FFP8 sim¬ 
ulations (Planck Collaboration XII 2016) using SMICAfor fore¬ 
ground separation. They are the same simulations we use to de¬ 


termine the error bars on our final SMICA results in Sect. 6. They 
contain the Collaboration’s best estimates of the CMB sky and 
of Planck's noise and beam effects, and have been cleaned by 
SMICA in the same way as the real sky map. The mask used is 
the same common mask defined for the real data analysis. For 
this test the estimators were all processed with the same settings 
used for the final data analysis. 

Here we take 159 of these maps, and process these using all 
the estimators. By contrast, for the final results in Sect. 6, the 
convergence of the error bars of each estimator was carefully 
checked, using more maps if required. This explains why there 
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Table 8. Results from the different estimators for /nl for the 
set of CMB simulations with = 12 described in Sect. 5.2. 
Results are given for T -only, £'-only, and the full combined T+E 
analysis. The shapes are assumed independent (see the main text 
for a discussion of this point). 




/nl 


Shape 

KSW 

Binned Modal 1 

Modal 2 


T 


Local . . . 

11.2 

+ 

6.7 

10.9 + 

6.3 

11.9 + 

6.6 

11.6 + 

6.6 

Equilateral 

26 

+ 

78 

24 + 

77 

31 

+ 

82 

27 + 

76 

Orthogonal 

-33 

+ 

34 

-33 + 

35 

-34 

+ 

36 

-33 + 

36 

E 

Local . . . 

11 

+ 

29 

12 + 

26 

9 

+ 

36 



Equilateral 

34 

+ 

182 

32 + 

153 

10 

+ 

241 



Orthogonal 

-37 

+ 

no 

-28 + 

115 

-31 

+ 

143 



T+E 

Local . . . 

11.3 

+ 

5.5 

11.2 + 

5.0 

11.1 

+ 

5.8 

11.0 + 

5.4 

Equilateral 

29 

+ 

52 

24 + 

50 

28 

+ 

54 

24 + 

50 

Orthogonal 

-29 

+ 

26 

-28 + 

23 

-30 

+ 

28 

-26 + 

23 


are some small differences between the error bars in Sect. 6 and 
the ones presented here. 

The results are shown in Table 9. Note that these are the re¬ 
sults from an independent analysis, without subtracting the ISW- 
lensing bias. We also show the results from Minkowski func¬ 
tionals (for the local case only).^^ We see that there is very good 
agreement between the bispectrum estimators even on this most 
complex and realistic set of simulations. The standard deviation 
of the difference between bispectrum estimators generally stays 
below one third of the error bar on /nl, the only exception being 
the T -only equilateral result for the Modal 1 pipeline, which is 
still smaller than one half of the error bar. We see that the results 
from Minkowski functionals are consistent, but clearly subopti- 
mal for /nl- They are however a valuable, independent check. 

The exact Fisher error bars for the nine shapes considered 
in the table are, in the same order as the table: 5.4, 69, 35; 31, 
131, 74; 4.7, 43, 21. Taking into account the relative error in the 
standard deviation of 1/ which is 5.6 % for 159 maps, 

we see that all bispectrum estimators are effectively optimal on 
all shapes, except for the £’-only equilateral case where they ap¬ 
pear slightly suboptimal. The small suboptimality of the Modal 
2 pipeline for the local shape seen here disappears once more 
maps are used (see the results in Sect. 6). 

In conclusion, all these validation tests show that we have 
very good agreement between the results from the different bis¬ 
pectrum estimators, not just on average, but also on a map-by- 
map basis. In addition we see that, despite the approximations 
made in the pipelines, and the simple treatment of the masked 
part of the maps (diffusive inpainting method and /sky factor), 
the bispectrum estimators we use are all essentially optimal. 


Since the Minkowski-functional pipeline automatically subtracts 
the ISW-lensing bias, the theoretical value for the bias as computed 
from the Fisher matrix has been added to its results, to make a direct 
comparison possible. 


6. Results 

6.1. Constraints on local, equilateral, and orthogonal fi^L 

In this section we investigate the local, equilateral, and orthogo¬ 
nal primordial templates. These are now established as the stan¬ 
dard shapes to study first when investigating the bispectrum 
(see Sect. 2 for a theoretical motivation and description of these 
shapes). However, they represent only the tip of the bispectral 
iceberg, and many more shapes are investigated in Sect. 8, while 
full model-independent reconstructions of the bispectrum are 
presented in Sect. 6.2. 

For a complete description of the Planck data set and the 
bispectrum estimator configurations we have used, we refer the 
reader in particular to Sect. 3.4 and Sect. 3.5. To summarize 
the overall analysis methodology, we have employed four inde¬ 
pendent bispectrum estimators on the full mission Planck tem¬ 
perature and polarization maps obtained from the four differ¬ 
ent component separation pipelines, SMICA, SEVEN, NILC and 
Commander. The bispectrum estimators are the KSW estima¬ 
tor with its skew-Q extension using exact separable templates 
(Sect. 3.1), the Binned estimator using fixed multipole bins 
(Sect. 3.3), and the Modal 1 and Modal 2 estimators, which both 
use separable eigenmode expansions (Sect. 3.2). Temperature is 
analysed over the multipole range ^min = 2 to ^max = 2000 or 
above and polarization is analysed from ^min = 40 to ^max = 
1500 or above (Sect. 3.5). By employing inpainting and a lin¬ 
ear term, all these estimators essentially achieve optimality (as 
shown by comparison with Fisher matrix forecasts). The linear 
term in Eq. (36) and the uncertainties are determined using the 
FFP8 simulations (Sect. 3.5), also processed through the dif¬ 
ferent foreground-separation pipelines. Our thorough validation 
campaign for these estimators is presented in Sect. 5. 

The results of the analysis of the four cleaned maps with the 
four estimators, for T-only, £’-only, and full T+E, are shown in 
Table 10, which is one of the main results of this paper. Results 
are determined while assuming all shapes to be independent, and 
are shown both with and without subtraction of the ISW-lensing 
bias (see Sect. 4.1 for more details about ISW-lensing). This bias 
is most important (relative to the size of the error bars) for the 
local shape, but also non-negligible for the orthogonal shape. 
Results here have not been marginalized over the point source 
contributions. While Sect. 4.2 shows that there is still a signifi¬ 
cant contamination by unclustered point sources in the cleaned 
maps, the correlation with the primordial templates is so small 
that this has no impact on the results reported here (as checked 
explicitly). 

While Table 10 is the main result of this Section, in order 
to simplify the use of the Planck results by the wider scientific 
community, we also present in Table 11 the results that can be 
considered the final Planck 2015 results for the local, equilateral, 
and orthogonal shapes. As in 2013, we select the combination of 
the KSW estimator and the SMICA map for this. The SMICA map 
consistently performs well in all data validation tests, which are 
discussed in detail in Sect. 7. The KSW estimator, while unable 
to deal with non-separable templates, treats separable templates 
exactly, and the local, equilateral, and orthogonal template are all 
separable. On the other hand, the binned and modal estimators 
can deal with non-separable shapes and have other advantages 
as well (like full bispectrum reconstruction), but at the price of 
using approximations for the templates. However, they have all 
been optimized in such a way that the correlation with the ex¬ 
act templates for the three primordial shapes is close to perfect, 
so that in the end the results by the different estimators are sta¬ 
tistically equivalent. Compared to the corresponding values in 
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Table 9. Results from the different estimators for /nl for the set of SMICA simulations based on FFP8 described in Sect. 5.3. Both 
the results for the estimators individually and for the differences with KSW are given, for T -only, £’-only, and the full combined 
T+E analysis. The shapes are assumed independent and the lensing-ISW bias has not been subtracted. 


/nl 

Shape KSW Binned Modal 1 Modal 2 Mink.F. B - KSW Ml - KSW M2 - KSW MF - KSW 

T 


Local. 

7.1 + 

5.5 

7.0 

+ 

5.4 

6.2 

+ 

5.5 

6.3 + 

6.2 

Equilateral . 

2 + 

67 

4 

+ 

67 

-4 

+ 

73 

5 + 

66 

Orthogonal . 

-23 + 

32 

-24 

+ 

33 

-24 

+ 

33 

-20 + 

36 

E 

Local. 

0.5 + 

32 

0 

+ 

35 

1 

+ 

30 



Equilateral . 

7 + 

144 

7 

+ 

143 

9 

+ 

152 



Orthogonal . 

5 + 

72 

7 

+ 

75 

4 

+ 

73 



T+E 

Local. 

5.6 + 

5.1 

5.0 

+ 

4.9 

4.7 

+ 

4.8 

4.3 + 

5.3 

Equilateral . 

3 + 

46 

5 

+ 

44 

3 

+ 

46 

4 + 

43 

Orthogonal . 

-10 + 

22 

-9 

+ 

22 

-9 

+ 

21 

-7 + 

22 


7 + 12 

-0.1 + 1.1 

-0.9 + 

1.9 

-0.8+ 1.9 

0 + 11 


2+19 

-6 + 

32 

3+18 



-1 + 11 

-0.9 + 

9.1 

3 +14 


0 + 49 

-0.8+ 8.3 

0.7 + 

8.3 


0 + 37 


0+37 

2 + 

35 




2 +22 

-1 + 

17 



5 + 11 

-0.6+ 1.2 

-0.9 + 

1.5 

-1.3+ 1.7 

-1 + 11 


2+14 

0 + 

14 

1.0+ 9.7 



0.8+ 7.0 

0.8 + 

7.3 

2.7+ 7.7 



Table 10, the difference in the numbers in the last column of 
Table 11 is due to the fact that in the latter the equilateral and 
orthogonal /nl have been determined jointly. 

Focusing on the results for temperature-only and the full 
temperature plus polarization {T+E) results, we see that there is 
no evidence for a non-zero bispectrum with any of these three 
primordial shapes (local, equilateral, and orthogonal). After 
ISW-lensing subtraction, all /nl values are consistent with 0 
at 68% CL. The temperature results are all very similar to 
the ones from the nominal mission data published in 2013 
(Planck Collaboration XXIV 2014), with very minor improve¬ 
ments in the error bars due to the additional temperature data. We 
also see that results are quite consistent when including polariza¬ 
tion, with error bars shrinking by about 15 % for local, 35 % for 
equilateral, and 40 % for orthogonal. 

Table 10 displays very good agreement between the results 
from the different estimators, at the level expected from the val¬ 
idation tests in Sect. 5. We also note how the error bars, which 
are determined using the FFP8 simulations, are statistically in¬ 
distinguishable from the optimal Fisher expectation. 

Different component separation methods also show a good 
level of agreement when looking at temperature-only and com¬ 
bined temperature plus polarization results. The accuracy of this 
statement will be shown and quantified in detail in Sect. 7. 
However, in the same section, we will also show how the agree¬ 
ment between /nl extracted from different cleaned maps be¬ 
comes significantly degraded when considering polarization- 
only resultsThe reasons behind this loss of internal consis¬ 
tency are not fully understood at present. Polarization data are, 
however, much noisier than temperature data, implying that the 
EEE bispectra have a close to negligible weight in the final 
combined measurement, which is dominated by the TTT and 
TTE configurations. In fact, as just mentioned above, the com¬ 
bined measurement looks perfectly self-consistent: local, equi¬ 
lateral and orthogonal /nl measurements in the T+E column of 
Table 10 pass all our tests of robustness. 

We can thus conclude that, while highly challenging from a 
technical point of view, the inclusion of polarization in our es- 

The F-only /nl agreement is still at a reasonable 1 cTf^^ level in 
most cases. However this is larger than expectations from simulations, 
as described in Sect. 7. 


timator pipelines has been a success, allowing for a significant 
tightening of the constraints on the three standard primordial bis¬ 
pectrum shapes. On the other hand, in light of the outstanding 
issues in E'-only results, we present our results conservatively, 
and recommend the reader to consider all /nl constraints that 
make use of polarization data throughout this paper as prelim¬ 
inary at the current stage. We stress again that this is a con¬ 
servative choice, which is made despite the fact that no test to 
date shows any evidence of leakage of the issues in EEE bispec¬ 
tra into the T+E measurements. A detailed description of all the 
data validation tests, which lead to the robustness-related assess¬ 
ments summarized here, can be found in Sect. 7 (for readers less 
interested in the technical details, the main results and conclu¬ 
sions of all these tests are summarized in Sect. 7.6). 

6.2. Bispectrum reconstruction 

6.2.1. Modal bispectrum reconstruction 

The starting point for modal bispectrum estimation is the ro¬ 
bust extraction of the modal coefficients jSn from each of the 
full mission foreground-separated maps, that is, SMICA, SEVEM, 
NILC, and Commander. The-coefficients are obtained for each 
of the temperature, polarization, and mixed bispectrum compo¬ 
nents, TTT, TTE, TEE, and EEE. Their cross-correlation between 
cleaning methods is an important validation of their accuracy, as 
we shall discuss in the next section, with excellent correspon¬ 
dence for temperature and some differences remaining in polar¬ 
ization. The modal basis number ^max = 2001 for the full mis¬ 
sion analysis has been substantially increased offering a higher 
effective resolution when compared to the 2013 Planck Data 
Release where ^max = 601 modes were used. Several different 
basis functions have been used, including trigonometric func¬ 
tions, sinlog basis functions, and polynomials (closely related to 
Legendre functions), with the latter chosen because of excellent 
convergence in the squeezed and flattened limits. 

We can reconstruct the full three-dimensional Planck bispec¬ 
trum, obtained using these basis functions, to visualize its main 
properties and to determine robustness. A comparison between 
the temperature-only bispectra from the Nominal Mission and 
full mission at the same ^max = 601 modal resolution is shown 
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Table 10. Results for the /nl parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW, binned 
and modal estimators from the SMICA, NILC, SEVEN, and Commander foreground-cleaned maps. Results have been determined using 
an independent single-shape analysis and are reported both without and with subtraction of the ISW-lensing bias; error bars are 68 % 
CL. 


/nl 


Independent ISW-lensing subtracted 

Shape KSW Binned Modal 1 Modal 2 KSW Binned Modal 1 Modal 2 


SMICA T 

Local. 10.2+ 5.7 8.7 + 5.4 6.8 + 5.5 7.8 + 6.0 2.5 + 5.7 1.3 + 5.4 0.5 + 5.5 1.7 + 6.0 


Equilateral. . 


-13 + 

70 

-26 + 

66 

-16 + 

67 

Orthogonal . 


-56 + 

33 

-41 + 

33 

-47 + 

33 

SMICA E 

Local. 


26 + 

32 

35 + 

34 

20 + 

30 

Equilateral . 


144 + 

141 

156 + 

143 

147 + 

159 

Orthogonal . 


-128 + 

72 - 

-128 + 

75 - 

-137 + 

73 

SMICA T+E 

Local. 


6.5 + 

5.0 

5.8 + 

4.9 

4.0 + 

4.8 

Equilateral . 


3 + 

43 

12 + 

44 

5 + 

48 

Orthogonal . 


-36 + 

21 

-34 + 

22 

-30 + 

21 

SEVEM T 

Local. 


11.3 + 

5.7 

9.7 + 

5.4 

8.1 + 

5.8 

Equilateral . 


—3 + 

69 

-16 + 

66 

-11 + 

75 

Orthogonal . 


-59 + 

33 

-47 + 

33 

-49 + 

34 

SEVEM E 

Local. 


60 + 

42 

62 + 

42 

44 + 

38 

Equilateral . 


292 + 

167 

320 + 

154 

302 + 

183 

Orthogonal . 


-184 + 

91 - 

-156 + 

93 - 

-172 + 

91 

SEVEM T+E 

Local. 


9.3 + 

5.2 

8.3 + 

4.9 

6.4 + 

5.0 

Equilateral . 


9 + 

47 

21 + 

48 

15 + 

52 

Orthogonal . 


-50 + 

23 

-46 + 

23 

-44 + 

23 

NILC T 

Local. 


10.5 + 

5.6 

8.7 + 

5.4 

6.4 + 

5.6 

Equilateral . 


-28 + 

69 

-45 + 

66 

—31 + 

75 

Orthogonal . 


-67 + 

33 

-48 + 

33 

-50 + 

33 

NILC E 

Local. 


0 + 

33 

18 + 

36 

-1 ± 

30 

Equilateral . 


75 + 

140 

97 + 

141 

64 + 

162 

Orthogonal . 


-79 + 

76 

-96 + 

81 

-78 + 

77 

NILC T+E 

Local. 


6.9 + 

5.1 

6.1 + 

4.9 

3.3 + 

4.9 

Equilateral . 


-9 + 

44 

-4 + 

44 

-15 + 

50 

Orthogonal . 


-35 + 

21 

-31 + 

22 

-27 + 

23 

Commander 

Local. 

T 

9.6 + 

6.1 

9.4 + 

5.7 

6.4 + 

6.6 

Equilateral . 


-19 + 

71 

-36 + 

68 

—3 + 

77 

Orthogonal . 


-49 + 

35 

-38 + 

34 

-49 + 

36 

Commander E 

Local. 

33 + 

39 

56 + 

40 

28 + 

37 

Equilateral . 


327 + 

165 

369 + 

157 

278 + 

178 

Orthogonal . 


-52 + 

88 

-70 + 

88 

-56 + 

81 

Commander 

Local. 

T+E 

7.7 + 

5.2 

7.9 + 

5.0 

5.2 + 

5.4 

Equilateral . 


16 + 

46 

26 + 

45 

30 + 

50 

Orthogonal . 


-37 + 

22 

-37 + 

23 

-39 + 

23 


-12 + 

68 

-11 

+ 

70 

-27 

+ 

66 

-12 

+ 

67 

-13 + 

68 

-63 + 

36 

-34 

+ 

33 

-14 

+ 

33 

-20 

+ 

33 

-44 + 

36 



26 

+ 

32 

34 

+ 

34 

20 

+ 

30 





144 

+ 

141 

155 

+ 

143 

147 

+ 

159 





-128 

+ 

72 - 

-126 

+ 

75 - 

-137 

+ 

73 



4.8 + 

4.9 

0.8 + 

5.0 

0.7 + 

4.9 

-0.6 + 

4.8 

0.7 + 

4.9 

6 + 

42 

3 

+ 

43 

9 

+ 

44 

3 

+ 

48 

5 + 

42 

-37 + 

21 

-25 

+ 

21 

-24 

+ 

22 

-21 

+ 

21 

-30 + 

21 

9.3 + 

6.0 

3.6 + 

5.7 

2.3 + 

5.4 

1.4 + 

5.8 

3.1 + 

6.0 

-6 + 

68 

-2 

+ 

69 

-18 

+ 

66 

-12 

+ 

75 

-7 + 

68 

-66 + 

36 

-36 

+ 

33 

-20 

+ 

33 

-23 

+ 

34 

-48 + 

36 



60 

+ 

42 

61 

+ 

42 

44 

+ 

38 





292 

+ 

167 

318 

+ 

154 

302 

+ 

183 





-183 

+ 

91 - 

-154 

+ 

93 - 

-172 

+ 

91 



7.9 + 

5.0 

3.3 + 

5.2 

2.8 + 

4.9 

2.1 + 

5.0 

3.5 + 

5.0 

5 + 

45 

8 

+ 

47 

17 

+ 

48 

14 

+ 

52 

4 + 

45 

-55 + 

22 

-39 

+ 

23 

-35 

+ 

23 

-33 

+ 

23 

-47 + 

22 

8.0 + 

6.2 

3.0 + 

5.6 

1.4 + 

5.4 

0.3 + 

5.6 

2.2 + 

6.2 

-15 + 

66 

-28 

+ 

69 

-47 

+ 

66 

-30 

+ 

75 

-17 + 

67 

-63 + 

35 

-45 

+ 

33 

-22 

+ 

33 

-28 

+ 

33 

-44 + 

35 



-1 

+ 

33 

17 

+ 

36 

-2 

+ 

30 





75 

+ 

140 

96 

+ 

141 

64 

+ 

162 





-78 

+ 

76 

-94 

+ 

81 

-78 

+ 

77 



5.3 + 

5.2 

1.2 + 

5.1 

0.9 + 

4.9 

-2.4 + 

4.9 

4.4 + 

5.2 

8 + 

42 

-9 

+ 

44 

-7 

+ 

44 

-16 

+ 

50 

4 + 

42 

-32 + 

21 

-25 

+ 

21 

-21 

+ 

22 

-16 

+ 

23 

—26 + 

21 

7.9 + 

6.3 

4.0 + 

6.1 

2.4 + 

5.7 

1.4 + 

6.6 

3.3 + 

6.3 

-14 + 

70 

-20 

+ 

71 

-38 

+ 

68 

-4 

+ 

77 

-18 + 

70 

-45 + 

37 

-29 

+ 

35 

-12 

+ 

34 

-25 

+ 

38 

-28 + 

37 



33 

+ 

39 

55 

+ 

40 

28 

+ 

37 





327 

+ 

165 

368 

+ 

157 

278 

+ 

178 





-52 

+ 

88 

-67 

+ 

88 

-56 

+ 

81 



6.8 + 

5.2 

3.7 + 

5.2 

3.0 + 

5.0 

1.6 + 

5.4 

3.7 + 

5.2 

29 + 

46 

14 

+ 

46 

23 

+ 

45 

28 

+ 

50 

26 + 

46 

-35 + 

22 

-29 

+ 

22 

-27 

+ 

23 

-30 

+ 

23 

-28 + 

22 


in Fig. 4. Note the excellent agreement with all the main fea¬ 
tures replicated in the new data. In Fig. 4 in the third bispec¬ 
trum, we also demonstrate the much higher bispectrum resolu¬ 
tion achieved with the full ^max = 2001 modes. The tetrapyd 
shape reflects the constraints on the wavenumbers T’l,^’ 2 , and T’s, 
with the squeezed configuration appearing on the axes that lie 


along one / = 0. The expected ISW-lensing bispectrum is an os¬ 
cillating signal in the squeezed limit along the tetrapyd edges; it 
is now measured with a significance of 3.0 cr (see Sect. 4.1). This 
ISW-lensing signal sets an interesting benchmark or threshold 
against which to compare the other strong features observed in 
the bispectrum and now defined with greater precision. The orig- 
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Fig. 4. Modal bispectrum reconstruction for Planck 2013 (top left) and 2015 (top right) temperature-only data, both using the SMICA 
maps. Here, we restrict the 2015 resolution to the same as 2013, using similar polynomials with ^max = 601. The two bispectra are 
very close to being in complete agreement in the signal-dominated regime shown up to ^max = 1500. In the lower panel, we show 
the Planck 2015temperature bispectrum at high resolution using the full ^max = 2001 polynomial modes. Large-scale features in 
the top panels become subdivided but the main 2013 signals remain, notably a stronger measurement of the ISW-lensing signal (the 
regular oscillations in the squeezed limit). 


inal “plus-minus” feature, with a large positive red peak around 
^ 150 followed by a larger negative peak near i ^ 250, re¬ 

mains though with more substructure, together with a broad neg¬ 
ative peak in the equilateral limit around I ^ 900, which can be 
associated with the third acoustic peak from the transfer func¬ 
tions. Oscillatory models, which can connect these three peaks. 


achieve higher significance. The apparent signal observed in the 
fiattened limit remains, with a distinct pattern of blue and red 
features on the surface of the tetrapyd. 

We also include a comparison with WMAP-9 in Fig. 5, 
where we have restricted the reconstructions to ^max = 600 for 
comparison with ^max = 601 modes. These plots, using identical 
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Fig. 5. Modal reconstruction for the WMAP-9 bispectrum (left) and the Planck SMICA 2015 T-only bispectrum (right) plotted for 
the domain i < 450, using identical isosurface levels. Here, we employed the full 2001 eigenmodes for both the Planck analysis 
at ^max = 2000 and for WMAP-9 analysis at ^max = 600, but for comparison purposes we have only used the first 600 eigenmodes 
in order to obtain a comparable resolution. The main features in the WMAP-9 bispectrum have counterparts in the Planck version, 
revealing an oscillatory pattern in the central region, as well as features on the tetrapyd surface. The WMAP-9 bispectrum has a 
much larger noise signal beyond ^ = 350 than the more sensitive Planck experiment, leading to apparent residuals in this region. 


Table 11. Results for the /nl parameters of the primordial lo¬ 
cal, equilateral, and orthogonal shapes, determined by the KSW 
estimator from the SMICA foreground-cleaned map. Both inde¬ 
pendent single-shape results and results with the ISW-lensing 
bias subtracted are reported; error bars are 68 % CL. The differ¬ 
ence between the last column in this table and the correspond¬ 
ing values in the previous table is that in the second column here 
the equilateral and orthogonal shapes have been analysed jointly. 
The final reported results of the paper are shown in bold. 


/nl(KSW) 


Shape and method Independent ISW-lensing subtracted 


SMICA (T) 

Local. 10.2 + 5.7 2.5 + 5.7 

Equilateral. -13 +70 -16 + 70 

Orthogonal. -56 + 33 -34 + 33 

SMICA {T+E) 

Local. 6.5 + 5.0 0.8 + 5.0 

Equilateral. 3 +43 -4 +43 

Orthogonal. -36 + 21 -26 + 21 


isosurfaces, show the same bispectrum structure including the 
“plus-minus” feature clearly bisecting the main i = 200 peak 
and the first oscillation of the ISW-lensing bispectrum visible 
along the lower tetrapyd edges. The WMAP-9 reconstruction 
only shows significant differences from Planck in the top right 
region, where the higher noise levels in WMAP-9 make its re¬ 
construction less reliable. 

All four components of the temperature and polarization 
bispectrum reconstruction obtained from SMICA are shown in 
Fig. 6. A direct comparison of the EEE polarization bispectrum 
for SEVEM, NILC and Commander, is shown in Fig. 7, where we 


note that these are orthogonalized E’-mode contributions (see the 
Modal 2 discussion in Sect. 3). It is interesting to observe pat¬ 
terns of features evident in the polarization bispectra from the 
different foreground-cleaned maps, which, although inherently 
noisier, have qualitative similarities. At a quantitative level, how¬ 
ever, the polarization bispectra modes from different methods are 
less correlated in polarization than in temperature, as we discuss 
in Sect. 7. 


6.2.2. Binned bispectrum reconstruction 

The (reconstructed) binned bispectrum of a given map is a 
natural product of the binned bispectrum estimator code (see 
Sect. 3.3). To test if any bin has a significant NG signal, we 
study the binned bispectrum divided by its expected standard 
deviation, a quantity for which we will use the symbol 
With the binning used in the estimator, the pixels are dominated 
by noise. We thus smooth in three dimensions with a Gaussian 
kernel of a certain width (Tbin- To avoid edge effects, due to the 
sharp boundaries of the domain of definition of the bispectrum, 
we renormalize the smoothed bispectrum, so that the pixel val¬ 
ues would be normal-distributed for a Gaussian map. 

In Figs. 8 and 9, we show slices of this smoothed binned 
signal-to-noise bispectrum with a Gaussian smoothing of 
(Tbin = 2, as a function of ii and ^ 2 - Very red or very blue regions 
correspond to a significant NG of any type. The two figures only 
differ in the value chosen for the 7 ’ 3 -bin, which is [518,548] for 
the first figure, and [1291,1345] for the second. We have de¬ 
fined two cross-bispectra here: + B^^^ -r B ^^^; 

and B^^^ = B^^^ + B^^^ + B^^^. These two cross-bispectra 
are then divided by their respective standard deviations (taking 
into account the covariance terms) to produce the correspond¬ 
ing and . Those three different permutations are not 
equal a priori due to the condition ii < 12 < h, which is imple¬ 
mented in the code to reduce computations by a factor of six. 


22 



















Planck Collaboration: Planck 2015 Results. Constraints on primordial NG 



Fig. 6. CMB temperature and polarization bispectrum reconstructions for Planck SMICA maps using the full set of polynomial modes with 
^max = 2001 and with signal-to-noise weighting. The top bispectra are the symmetric pure temperature TTT (left) plotted with £ < 1500 and 
£^-mode polarization FEE (right) shown for 30 < ^ < 1100. Below are the mixed temperature/polarization bispectra with TTE on the left (with E 
multipoles in the z-direction) and TEE on the right (with T multipoles in the z-direction). All S/N thresholds are the same. 



Fig. 7. Comparison of CMB polarization bispectrum FEE reconstructions for Planck NILC, SEVEM, and Commander foreground-separated maps 
with signal-to-noise weighting. Note that these results are not as internally consistent between the four methods, also comparing SMICA shown in 
Fig. 6, which is closest to NILC. We will compare the underlying modal coefficients below to demonstrate these differences quantitatively. 
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Fig. 8 . Smoothed binned signal-to-noise bispectrum B for the Planck 2015 cleaned sky map, as determined with the binned estima¬ 
tor, as a function of and I 2 for a fixed ^ 3 -bin [518, 548]. From left to right results are shown for the four component separation 
methods SMICA, SEVEM, NILC, and Commander. From top to bottom are shown: TTT, TTT cleaned from radio and CIB point sources; 
T2E, TE2; and EEE. The colour range is in signal-to-noise from -4 to -1-4. The light grey regions are where the bispectrum is not 
defined, either because it is outside the triangle inequality or because of the cut = 2000 . 


However, part of the smoothing procedure involves adding the 
other five identical copies, so that in the end the plots are sym¬ 
metric under interchange of and ^2 (and is symmetric 
under interchange of all its indices). The grey areas in the plots 
are regions where the bispectrum is not defined, either because 


it is outside of the triangle inequality, or because of the limita¬ 
tion - 2000. Given that in both plots ^3 is fixed at less than 
2000 , is not defined if both ii and I 2 are larger than 2000 , 
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while is undefined if either or £2 (or both) are larger 
than 2000 . 

Results are shown for the four component separation meth¬ 
ods SMICA, SEVEN, NILC, and Commander, and for TTT, T2E, 
TE2, and EEE. In addition we show on the second line of each 
figure the result for TTT with the radio (unclustered) and CIB 
(clustered) point source bispectra subtracted according to their 
jointly measured amplitudes. It is clear, in particular in the sec¬ 
ond figure, that at higher i there is a very significant point source 
contamination in the cleaned TTT bispectra, in agreement with 
the results of Table 3. However, after removing it we do not see 
a clear signal of any other residual NG. Of course this is for the 
moment only a qualitative statement; more quantitative tools for 
studying the amount of NG in these smoothed bispectra are in 
development. 

Looking at the polarized bispectra in the high -^3 slices, in 
particular TE2 and EEE, we do see some bluer and redder re¬ 
gions that might indicate residual NG. This agrees with state¬ 
ments made earlier, and discussed in greater detail in the next 
section, that the Planck polarized bispectrum is for the moment 
not as clean and well-understood as the temperature one. We also 
see a very good qualitative agreement between the four compo¬ 
nent separation methods in temperature, which worsens some¬ 
what when mixing in more and more polarization; in particular 
SMICA and NILC give very similar results. 

6.3. Primordial Curvature Reconstruction 

In this section, we compress the information in T and E maps 
into maps of projected primordial curvature fiuctuations, f. To 
the extent that the primary CMB temperature and polarization 
are Gaussian fields and the early universe f fiuctuations are pre¬ 
sumed Gaussian, all of the information in T, E is encoded in 
the mean field (Wiener-filtered map) plus fiuctuations, character¬ 
ized by a covariance matrix; it is just a different compression/re¬ 
expression of the original statistics. Thus, one could look for 
non-Gaussian deviations by directly evaluating the 3-point and 
higher statistics of this re-expression. Since fiuctuations as well 
as the mean are given, a full description of the errors follows. 
That a Gaussian assumption is made does not mean that the 
Wiener-filtered map is in fact Gaussian, since the constraining 
data may drive it from a map consistent with Gaussianity. Hence 
one can search in the statistics of the Wiener-filtered map for 
evidence for non-Gaussianity. 

The weighting of the temperature and £’-mode polarization 
associated with the Wiener filter is an optimal one (inverse- 
total-covariance weighted), and estimators for non-Gaussianity 
constructed using this expression are nearly optimal. However, 
since a (fiducial) primordial power spectrum is assumed for the 
Wiener filter (usually a uniform Us ^-spectrum, using the best fit 
Planck 2015 scalar spectral index and power amplitude As), this 
approach is best suited for perturbative non-Gaussianity of the 
sort we treat in this paper. In practice, the estimators used here 
for non-Gaussianity act directly on temperature and polarization 
data, rather than through the intermediary of optimal f-maps and 
their fiuctuations. 

The scalar fiuctuations can be expressed in terms of the cur¬ 
vature variable f(x) = lna(x) on uniform total density hypersur¬ 
faces, where a is the inhomogeneous expansion factor, or, equiv¬ 
alently, by its wavenumber transform f(k). In turn, f(k) can be 
expanded in multipoles, ^LuAk). Instead of the magnitude of the 
wavenumber k, a mixed representation gives a multipole expan¬ 
sion at each comoving distance from our location, ^lmk(x)^ 


M = 0,..., L. Here k = c or s, with k = c referring to the real part 
of ^LM and K = s io the imaginary part. For M = 0 there is no 
s component, only c. The mean of the f field given the temper¬ 
ature field T and its covariance describing allowed fiuctuations 
about that mean are 

I «Im.> = cf Or) [cy ]■’ ( 52 ) 

{^^LMKix)^^L'M'K'ix')\T) = ^LL'^MM'^kk' X (53) 

{cf(Y,V) - cf{x)[cVTcl^^')]. 


In these expressions, is the total covariance for the tem¬ 
perature, including both signal and noise variances. The specific 
forms assume the matrix is diagonal in multipole space, 
depending only on , but if the noise is inhomogeneous it 
will have off-diagonal components and the equations become 
matrix equations. In this section, we assume the noise is diag¬ 
onal. Eq. (52) shows the unsurprising result that for each LMk 
only one mode is determined by T. Replacing T hy E gives 
{^lmk(x) I Although T and E are correlated, 

the uncorrelated part of E delivers a different mode for f from 
the one given by T. Thus when f is constrained by both T and 
E, the two modes deliver substantially more information than for 
T and E alone, and the fiuctuations about the mean are thereby 
diminished. This is further helped by the acoustic oscillations of 
polarization being out of phase with those for T, so when T is 
down, E is not, and vice versa. The interplay of the two modes 
is quantified by signal-to-noise in Fig. 10. 

When both T and E are included as constraints, a two by 
two matrix appears which includes the correlation in the 
off-diagonal as well as the and along the diagonal: 


{^LMk(x) I ^LMk^ ^LMk) - 

. rrTT 

[cf(Y) cf(Y)l 




TE 

cr ci^ 


-1 r T 
^LMk 


^LMk 


\T,E) - SLLrduM'^KK'iC^iiX^X') 




r^TT r^TE 
f^ET /^EE 


-1 


cfOr') 
cf Or') J 


(54) 


(55) 


Fig. 11 shows an all-sky mean field (Wiener-filtered) map of 
the curvature variable f constructed in (densely-packed) shells 
from the multipoles ^lmk(x)- The figure actually shows f pro¬ 
jected onto differential visibility, w = de~'^ Idx- The map fy;(^, 0) 
is the spherical transform of ^lmkw = f ^(x)^X^lmk(x)’ The top 
left map is {^u)\T) using the all-sky maps SMICA DXl 1. The mid¬ 
dle left map is {^uj\E), based on just the polarization information. 
(The E-maps did not have high pass filtering imposed, so that 
an indication of the full results obtainable with E alone can be 
seen.) The bottom left panel shows E). Visually the com¬ 
bined two-mode constraint has more detail than either of the one¬ 
mode constrained maps, enhancing the estimation of fy;. In each 
of the three cases, two sample realizations are shown which in¬ 
clude allowed fiuctuations about the mean map. These illustrate 
the level of uncertainty in the patterns found in the mean maps, 
necessary ingredients to the full statistical description. Note that 
the fiuctuations about the mean are larger for T or E alone than 
for the combined T, E. 

In all cases, the fiuctuations become very large indeed if we 
allow for arbitrary projections, not conditioned by the visibility. 
Another approach similar to the differential visibility projection 
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Fig. 9. Similar to Fig. 8, but with €3 e [1291,1345]. 


is to make a slice at;^ = Xdec'^ making the projection a delta- 
function in;^ at the peak where dw/d Inx = 0, which defines 
XdQc- This single slice case is shown in Fig. 12. The figures look 
quite similar in structure to the projected differential visibility 
plots. Although;^ ~ Xd&c is where most information is, the fluc¬ 
tuations of a single slice are somewhat larger than for the differ¬ 
ential visibility projection. The differential visibility includes a 
subdominant “reionization bump”, with a small weight focused 
at low L where sample variance is large; restricting the projec¬ 


tion of f to be just over the reionization region results in a low 
L mean map that is largely swamped by the allowed fluctuations 
about it. The fluctuations about the mean become much larger if 
we project over all;^, since there are vast terrains in;^ in which 
there is very little CMB temperature or polarization information; 
in particular between recombination and reionization. For that 
uncharted territory, a realization of the reconstructed f-map re¬ 
verts to a realization of the fiducial f-power. 
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Temperature S/N 




X/Xre 


components off-diagonal in;^): 


{smiix) 


P = 


ci^(x) 




cvcf{x,x) 


(57) 


for r, with an equivalent expression for E. The signal-to-noise 
structure in can be contrasted with the differential visibility 
plotted in the bottom panel: Fig. 10 shows that the L-x in¬ 
formation has a reach beyond the differential visibility struc¬ 
ture, especially for low multipoles, because the associated waves 
can straddle the last scattering surface. The “reionization bump” 
in signal-to-noise seen in Fig. 10 is not that prominent. The 
Integrated Sachs Wolfe impact on the S/N is evident, but is rel¬ 
atively low, with the consequence that we cannot draw out high 
significance results from the ISW effect for f reconstruction (or 
equivalently for gravitational potential reconstruction). 

The top panels in Fig. 10 are for an ideal experiment, with 
no noise (apart from the cosmic variance “noise” in the Wiener 
map fluctuations). The bottom panels of Fig. 10 include realistic 
Planck noise, as estimated from the FFP8 simulations: the high 
signal to noise in E in the top panels is noticeably diminished 
around L ~ 100 over what a cosmic variance limited experiment 
would give. 

Figs. 11 and 12 are Gaussian-filtered on a relatively large 
L ~ 200 scale (40' FWHM). To see what happens at higher 
resolution. Figs. 13 and 14 zoom in on a typical 20° x 20° patch, 
with long waves removed using a filter Wl = sin^(L-Lc)/AL for 
Lc < L < Lc-\-AL. The specific choices are = 20 and AL = 20, 
the values used for the Planck high pass for polarization maps. 
To allow direct comparison, we have just done the same high 
pass filtering for the T map. We also removed the means of the 
maps. The resolution is L ~ 400 (20' FWHM). Figs. 13 and 
14 illustrate that the fluctuations play a larger role in this higher 
resolution regime. Note that the Xdtc slice has about the same 
fluctuation level as the differential visibility projection. 

Tensing effects are not taken into account in these f-maps. 
In principle one could de-lens the temperature and polarization 
maps before forming the Wiener filter. In practice this would be 
a highly noisy operation with current Planck data. Hence the 
contaminating influence of lensing on these f-reconstructions 
would be treated by comparing simulated f-maps with and with¬ 
out lensing. Such corrections are expected to be a subdominant 
bias, as it has been modelled in the rest of the paper. 


Fig. 10. Signal-to-noise for T and as a function of multipole 
L and comoving distance for an ideal cosmic-variance lim¬ 
ited experiment (upper panels) and for the Planck experiment, 
with noise determined from FFP8 simulations (lower panels). 
The complementary nature of the information provided by T and 
E is evident. The differential visibility is shown for comparison. 


Fig. 10 quantifies where the information resides in L-x 
space. Here, the signal variance is the ensemble average of the 
square of the mean field, assuming the f are drawn from a 
Gaussian with fiducial covariance, 

{{Umk(x) I T){^lmM') I T)) = C^,^{x)[CYr'cl^(x'\ (56) 

and the noise (i.e., fluctuation) variance is the covariance matrix, 
as given by Eq. (54) for T and the equivalent for E. Fig. 10 plots 
the signal-to-noise at each individual slice (hence ignores the 


7. Validation of Planck resuits 

In this section, we perform a battery of tests aimed at verify¬ 
ing the robustness of the results obtained in the previous sec¬ 
tion. Table 10 shows excellent agreement with our 2013 analy¬ 
sis of nominal mission data (Planck Collaboration XXIV 2014). 
The agreement using different component separation methods in 
temperature is also generally very good. Our focus here is thus 
on polarization bispectra. Redundancy is perhaps the most im¬ 
portant element in our analysis, as far as robustness is concerned. 
We devote considerable attention to comparing the outcomes of 
different estimators and component separation pipelines, and as¬ 
sess their level of internal consistency. We also verify the stabil¬ 
ity of our results in the harmonic and pixel domains, by consid¬ 
ering different sky cuts and multipole intervals. Given the large 
computational requirements of these tests, and since results from 
different optimal estimators agree very well, as shown in previ¬ 
ous sections, we will variously use the KSW, binned, or modal 
pipeline for different tests. In doing this we will also exploit 
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Fig. 11. Mean field and fiuctuations for differential visibility projected f, as described in the text. The filter used in these maps is a 
40' FWHM Gaussian. The grey regions are masked out. 


the complementarity of different decompositions, which might 
make some of them more suited for different tests than others 
(for example, the binned pipeline directly works with a harmonic 
space decomposition of the bispectrum, thus making it perfectly 
suited for tests of ^-dependence, the modal pipeline can perform 
quick model-independent tests by working on a relatively small 
subset of bispectrum modes, and so on). 

7.1. Dependence on foreground cleaning method 

7.1.1. Comparison between /nl measurements 

In Table 10 we show /nl results for the local, equilateral and 
orthogonal shapes, using four different optimal estimators, and 
four different foreground cleaning pipelines. The agreement be¬ 
tween different estimators, on a given map, is within a fraction 
of a standard deviation, in line with theoretical expectations and 
simulations, as reported in Sect. 5. This level of agreement ap¬ 
plies to all of the T+E, {TTT), and {EEE) bispectra. 

The overall picture becomes more complex when comparing 
outputs across different foreground cleaning methods and esti¬ 


mators. Whereas for TTT and T+E results the agreement also 
seems quite good in this case (being at the level of half a sigma or 
better, for nearly all combinations of cleaned maps and shapes), 
larger discrepancies are present in the EEE bispectrum measure¬ 
ments. The most notable differences are found for the equilateral 
shape, where SMICA and NILC find values of /nl consistent with 
0 within 1 cr, while SEVEN and Commander measure a roughly 
2 (T deviation from Gaussianity. The largest discrepancy is found 
for the pair Commander-NILC, using the binned pipeline (see ta¬ 
ble 10). This estimator recovers = 369 + 160 using the 

Commander £’-map and = 97 + 141 using the NILC £’-map. 
Other pipelines, and different choices of component separation 
methods, show slightly smaller but similar discrepancies, at a 
level of about l.Scr. The same shape and estimator analysis of 
temperature maps shows good agreement: Commander recovers 
= -36 ± 73, while NILC gives = -45 ± 71. The 
combined T+E measurement, still for the same modal pipeline 
and equilateral shape, yields = 26 + 50 for Commander and 
~ “4 + 46 for NILC, corresponding to about half a standard 
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Fig. 12. Mean field and fiuctuations for a Xdec slice of f, as described in the text. The filter used in these maps is a 40' FWHM 
Gaussian. The grey regions are masked out. 


deviation difference. This general trend is seen for other shapes 
and estimators. 

Simulations were used to give insight into the expected level 
of disagreement. For each of the four component separation 
methods, we generated 100 FFP8-based Gaussian simulations 
with realistic beams and noise levels. These simulations start 
from the same initial single frequency realizations, and are pro¬ 
cessed through the four different foreground cleaning pipelines. 
The starting maps do not include any foreground component (the 
same map generation procedure is used in the Monte Carlo deter¬ 
mination of error bars). The differences in final simulations are 
thus generated only by the different data filtering and coadding 
operations performed either in pixel, harmonic, or needlet do¬ 
mains by the various foreground cleaning methods, and by addi¬ 
tional manipulations of the maps which are required for /nl es¬ 
timation, such as inpainting. Therefore, the average scattering in 
/nl, measured from these realizations, provides us with a base¬ 
line assessment of the expected discrepancies between different 
methods when foreground residuals and other spurious sources 
of NG are negligible. We can then compare them with differ¬ 


ences observed on the data to establish whether they are consis¬ 
tent with expectations, or are too large. The latter would raise 
the concern that foreground contamination, or other systematics, 
might be affecting the results. 

Results are shown in Table 12, for EEE and T+E and two dif¬ 
ferent sky coverages. The scatter between /nl values from sim¬ 
ulations is about 0.5 (T for both T+E and EEE. This is smaller 
than the differences in the Planck /nl values obtained from EEE 
analysis of different foreground-cleaned maps, especially for the 
equilateral shape. However, for the final combined T+E mea¬ 
surement, observed differences are in good agreement with ex¬ 
pectations from simulations for the majority of cases. Another 
important point is that the consistency shown in Table 12 for 
T+E measurements is stable to the change of sky coverage (in 
polarization) from /sky = 0.74 to /sky = 0.64. This will be con¬ 
firmed by additional tests later in this section. For the SMICA- 
SEVEM pair we also verified stability using an even larger mask 
with /sky = 0.52. 

Residual foregrounds may be responsible for at least some 
of the observed excess of scatter in £’£’£’-derived /nl between 
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Fig. 13. Mean field and fiuctuations for differential visibility projected for a 20 degree by 20 degree patch, as described in the 
text. The filter used in these maps is a 20' FWHM Gaussian. Sources in grey are masked out. 


different cleaning algorithms. This is supported by the fact that 
several FEE results for this test change significantly for different 
masks, and that discrepancies are alleviated by using a larger 
mask, especially for equilateral shapes (see e.g., SEVEM-NILC 
and SMICA-SEVEM in Table 12). However, modal coefficients and 
their correlations are stable to a change of mask (see below), as 
are values of /nl for a given component separation method (see 
Table 13). 

Another possible contributor is a mismatch between the 
noise model (used to build the estimator normalization, weights 
and linear term), and the actual noise in the data. Polarization 
data are very noisy, and it is a known problem that the model 
assumed underestimates the true noise. This means that the error 
bars for FEE /nl results, quoted in Table 10, are somewhat un¬ 
derestimated, which does not seem to be a problem for the final 
T+E results, since the weight of the FEE bispectrum in the final 
combined measurement is very low. This is confirmed by the re¬ 
sults of this test. Indeed, we investigate FEE in detail because it 
is a useful and sensitive indicator of various systematics in the 
polarized maps (which could eventually leak into the TTE and 
TEE bispectra), rather than for its statistical weight in the final 
measurement. It is then fair to say that issues in the FEE bispec¬ 
tra, and related /nl measurements are not yet fully understood 
and will require further investigation. Even though the T+E are 


consistent, we recommend that results that include polarization 
data are regarded as preliminary at this stage. 


7.1.2. Comparison between reconstructed bispectra 

It is important to stress that the conclusions reached at the end 
of the previous subsection refer to the three main bispectra in 
our analysis, defined by the standard scale-invariant local, equi¬ 
lateral, and orthogonal primordial shapes. These shapes select a 
specific subset of configurations in the overall bispectrum do¬ 
main (essentially squeezed, equilateral, and flattened triangles). 
Therefore, testing consistency between methods for these shapes 
does not guarantee that results for the many other NG models 
considered in this work (such as, for example, the oscillatory 
bispectra of Sect. 8) will display the same level of agreement. 
For this reason we decided to perform a model-independent 
test of consistency between methods, based on comparisons be¬ 
tween the Pn eigenmodes used for bispectrum reconstruction 
in Sect. 6.2. We also reconstruct the bispectrum starting from 
a binned ^-decomposition, and this will be used in Sect. 7.4 
to study stability of the results in the harmonic domain. For 
the Pn study we consider a simple test based on measuring the 
correlation coefficient between modes extracted from different 
foreground-cleaned maps. The correlation is defined, as usual. 
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Fig. 14. Mean field and fiuctuations for a^dec slice of for a 20 degree by 20 degree patch, as described in the text. The filter used 
in these maps is a 20' FWHM Gaussian. Sources in grey are masked out. 


by 



cov 

((^n)i(0-n)j 


(58) 


and we measure it for each combination of the SMICA, SEVEN, 
NILC, and Commander maps, labelled by the indices ij. Results 
are given in Table 14 and 15 for the two modal pipelines and 
are illustrated in Fig. 15. These results show an excellent de¬ 
gree of correlation between different maps in temperature (es¬ 
pecially for SMICA, SEVEM, and NILC), which is reduced when 
polarization is considered. In fact the correlation for polarization 
is not much lower than temperature for SMICA and NILC, while 
it reduces the correlation for the pairs SMICA-SEVEM, and NILC- 
SEVEM, and for Commander when paired with any other method. 
This is consistent with previous findings of our /NL-based test. 
To test if these results are due to foreground residuals (or other 
effects that are not included in the simulations), we evaluate the 
same mode-mode correlations on the same sets of 100 realis¬ 
tic, foreground-free, Gaussian simulations as previously used, 
and processed through each of the different component separa¬ 
tion pipelines. For this analysis we consider TTT and EEE bis¬ 
pectra, expanded via the low-resolution modal estimator. Our 
results are reported in Table 14, in the simulation column, and 
they clearly show that the trend in the simulations is consistent 
with what we see in the Planck data. In particular, EEE results 


show a lower degree of correlation in simulated maps, for the 
same pairs of methods. The observed loss of correlation in po¬ 
larization does not seem to come from unresolved foregrounds 
or other unaccounted systematics, but rather something intrin¬ 
sic to the foreground-removal algorithms. They are substantially 
different, as SMICA and NILC both perform the cleaning in har¬ 
monic space, at the level of E and B multipoles, whereas SEVEM 
is essentially a pixel-space template fitting method, performing 
the subtraction on Q and U maps, or inpainted before /nl esti¬ 
mation. These issues will be studied in greater detail in future 
work, using Wiener-filtered, as well as inpainted maps for /nl 
estimation. However, we have already seen that the larger scat¬ 
ter between modes from different foreground cleaning methods 
does not have a serious impact on /nl estimation, at least for 
the standard local, equilateral, and orthogonal shapes. The non¬ 
standard shapes need to be analysed separately to check robust¬ 
ness of NG polarization results. This is the approach we will take 
for the various non-standard NG models. 


7.2. Dependence on sky coverage 

For each of the four component separation methods, we have 
used two different polarization masks, namely the same polar¬ 
ization mask as in Sect. 6, with /sky = 0.74 (defined as the po¬ 
larization “common mask” in Sect. 3.4), and an extended mask 
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Table 12. Comparison between local, equilateral, and orthogonal /nl results, obtained using the four different component separation 
pipelines. For each pair of cleaning methods, and for each NG model, we compute the difference in the measured /nl- The quoted 
error bar is the standard deviation of the same difference, extracted from a set of 100 realistic Gaussian simulations per method, not 
including foregrounds. These results have been obtained using the low resolution modal pipeline. See main text for comments and 
further details. 







/nl (methodi) 

- /nl (methoda) 








/sky 

= 0.74 



fsky 

= 0.64 



Methods 

Local 

Equilateral 

Orthogonal 

Local 

Equilateral 

Orthogonal 

SMICA-SEVEM 












T . 

-1.2+ 0.9 

-6.0 + 

8.7 

1.5+ 4.8 







E . 

-19 +21 

-155 

+ 

86 

34 +57 

5 +22 

-82 

+ 

90 

-11 

+ 66 

T+E . 

-2.4+ 1.6 

-10 

+ 

18 

13.5+ 9.4 

-1.5+ 1.7 

-12 

+ 

18 

13 

+ 10 

SMICA-NILC 












T . 

0.4+ 1.0 

14.5 + 

8.9 

2.5+ 4.7 







E . 

26 + 11 

83 

+ 

52 

-59 + 27 

26 + 13 

32 

+ 

56 

-96 

+ 28 

T+E . 

-0.7+ 0.9 

20.0 + 

8.2 

-3.3+ 3.8 

0.6+ 0.9 

18.4 + 

8.4 

-4.5+ 4.0 

SMICA-Commander 

T . 

0.4+ 3.5 

-14 

+ 

23 

1.7 + 14 







E . 

-3 + 16 

-130 

+ 

77 

-81 +42 

-13 + 17 

-117 

+ 

100 

-59 

+ 40 

T+E . 

-1.3+ 3.2 

-25 

+ 

18 

9+10 

-1.4+ 3.3 

-26 

+ 

18 

13 

+ 10 

SEVEM-NILC 












T . 

1.6+ 1.0 

20 

+ 

12 

1.0+ 4.5 







E . 

45 +26 

239 

+ 

94 

-94 + 69 

30 +29 

114 

+ 

105 

-86 

+ 79 

T+E . 

3.1+ 1.8 

30 

+ 

18 

-17 + 10 

2.2+ 1.9 

30 

+ 

18 

-18 

+ 10 

SEVEM-Commander 

T . 

1.6+ 3.4 

-8 

+ 

22 

0+14 







E . 

16 +22 

25 

+ 112 

-116 +59 

-18 +25 

-35 

+ 

121 

-48 

+ 64 

T+E . 

1.2+ 3.3 

-14 

+ 

21 

-5 + 11 

0.2+ 3.4 

-14 

+ 

20 

0 

+ 11 

NILC-Coininander 

T . 

0.0+ 3.0 

-28 

+ 

22 

-1 + 12 







E . 

-29 + 21 

-213 

+ 

84 

-22 + 54 

-39 + 23 

-149 

+ 

108 

38 

+ 55 

T+E . 

-1.9+ 3.1 

-45 

+ 

18 

12 + 11 

-2.0+ 3.2 

-44 

+ 

17 

18 

+ 11 


with /sky = 0.64. The temperature mask is kept unchanged in 
this test, and it covers a sky fraction /sky = 0.76 (the temper¬ 
ature “common mask” of Sect. 3.4). We report the variation in 
/nl for the three standard shapes in Table 13, which shows in¬ 
sensitivity to /sky, in agreement with earlier results on T+E. In 
this case, however, the EEE results also seem quite stable, sup¬ 
porting the view that foreground residuals are not affecting our 
local, equilateral, and orthogonal /nl results, especially for the 
final, combined T+E measurements. Tests on FFP8 simulations 
including foregrounds (see Sect. 7.3) suggest that /nl measure¬ 
ments obtained from the SMICA and SEVEN maps are the most 
accurate under the current choice of mask. As a further check of 
these two methods we consider a third polarization mask, with 
/g]^y = 0.53, and repeat the combined T+E /nl measurement, 
also finding stable results. For SMICA we find /^^^^ = 5.6 + 5.4, 
= 65 + 58, and = -30 ± 26, while for SEVEM we 
obtain = 9.4 + 5.4, “ = 75 + 59, ° = -50 + 30. 


We also perform model-independent checks by looking at 
the correlation coefficient between different sets of bispectrum 
modes, in a similar way to Sect. 7.1.2, but now changing the po¬ 
larization mask. Results are reported in Table 16, and confirm 
that the data and simulations behave similarly, and that polariza¬ 
tion modes display a lower correlation level than temperature. 


7.3. Tests on simulations 

We consider two realistic data simulations, one of which is 
Gaussian, while the other includes local NG. We start with a 
foreground-free realization, add foregrounds according to the 
Planck Sky Model, and finally process through the four com¬ 
ponent separation pipelines. By estimating /nl in the input 
foreground-free simulation, for each method, and comparing to 
/nl recovered from the cleaned maps (or with the input local /nl, 
for the NG case), we can assess both the impact of foregrounds 
on our measurement before subtraction and which method gives 
the highest accuracy. The necessity to clean is very apparent in 
the middle set of columns in Table 17, where no cleaning has 
been performed. 

SMICA and SEVEM give the best results, both in the G and NG 
tests. In the G test, reported in Table 17, SMICA results show an 
agreement between the input and the cleaned map at the level 
crf^^/2 for all shapes, and for all of TTT, EEE, and T+E. SEVEM 
displays a similar level of accuracy, except for where the 
difference is larger, but within one standard deviation. NILC and 
Commander clearly perform worse for the local shape, with NILC 
showing a 2crf^^ difference, and Commander even larger than 
that. In the NG test, reported in Table 18, SEVEM gives the most 
accurate results, recovering the input with o-f^^/2 accuracy or 
better. Results for SMICA are accurate at the 1 cTf^^ level for the 
TTT constraint, and worse (about 2 0 -/^^ ) EEE. However, the 
combined T+E measurement is again very good {(Tf^^l2). NILC 
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Table 13. For each of the four foreground cleaned maps, we compute /nl for the local, equilateral, and orthogonal modes using 
two different polarization masks, one with /sky = 0.74 and the other with /sky = 0.64, while for temperature we use a single mask 
with /sky = 0.76. We then calculate the difference between the two measurements and compare with expectations from simulations, 
obtained in the following way: firstly, we generate realistic Gaussian realizations for each component separation pipeline, not 
including foregrounds; then, for each simulated map and NG model, we measure /nl using the two masks in turn; finally, we 
calculate the standard deviation on 100 realizations. See the main text for more details and a discussion of these results, which were 
obtained using the low resolution modal pipeline. 




/sky = 0.74 



/sky = 0.64 



Difference 



Local 

Equilateral 

Orthogonal 

Local 

Equilateral 

Orthogonal 

Loeal 

Equilateral 

Orthogonal 

SMICA 

T . 

E . 

T+E . . . 

6.8+ 5.4 
25 +30 
4.0+ 4.8 

-17+ 66 
147 + 159 
5+46 

-48 + 33 
-137 + 73 
-30 + 21 

48 +31 
4.6+ 5.2 

220 + 168 
19+ 55 

-180 + 81 
-37 + 22 

-23 + 16 

-0.7+ 1.2 

-73 + 68 
-14 + 14 

43 +34 

6.7+ 7.7 

SEVEM 

T . 

E . 

T+E . . . 

8.1 + 5.8 
44 +38 
6.4+ 5.0 

-11 + 75 
302 + 183 
15+ 52 

-49 + 34 
-172 + 91 
-44 + 23 

43 +39 
6.2+ 5.3 

303 + 191 
31 + 54 

-170 + 96 
-50 + 25 

1 + 19 
0.2+ 1.3 

0 + 76 
-16 + 15 

-2 +49 
6.3+ 8.8 

NILC 

T . 

E . 

T+E . . . 

6.4+ 5.8 
-1 +30 

3.3+ 4.9 

-31 + 76 
64 + 162 
-15+ 50 

-50 + 33 
-78 + 77 
-27 + 23 

22 +30 
4.0+ 5.3 

190 + 162 

1 + 56 

-84 + 77 
-33 + 23 

-23 + 16 

-0.7+ 1.3 

-124 + 67 
-16 + 13 

6+37 
5.4+ 7.5 

Commander 

T . 

E . 

T+E . . . 

6.4+ 6.6 
28 +37 
5.2 + 5.4 

-3+ 77 
278 + 178 
30+ 50 

-49 + 36 
-56 + 81 
-39 + 23 

61 +38 

6.0+ 5.7 

337 + 188 
45+ 55 

-122 + 91 
-51 +25 

-32 + 20 
-0.7+ 1.5 

-60 + 92 
-15 + 14 

66 +47 
11.5+ 8.9 


Table 14. Correlation coefficients between pairs of bispectrum modes, extracted using two different component-separated maps. For 
both TTT and EEE we compare correlations measured from data with averages over 100 Gaussian realizations. The simulations are 
processed through the different component separation pipelines in the same way as the data, but they do not include any foregrounds. 
The correlation is clearly lower for EEE bispectra than for TTT. However, this is seen not only in data but also in simulations, 
indicating that it is not due to foreground residual contamination or other unaccounted for systematics. The results presented in this 
table are obtained using the low resolution modal pipeline, with 610 modes; results on data have also been cross-checked with the 
high-resolution modal estimator, using 2001 modes, and they are stable. 


Methods 


/sky — 

0.74 


/sky 

= 0.64 

77T 

Data 

Simul. 

EEE 

Data 

Simul. 

EEE 

Data Simul. 

SMICA-SEVEM . 

0.97 

0.97 

0.61 

0.62 

0.60 

0.61 

SMICA-NILC . 

0.97 

0.97 

0.95 

0.95 

0.95 

0.95 

SMICA-Commander ... 

0.78 

0.81 

0.70 

0.70 

0.73 

0.73 

SEVEM-NILC . 

0.96 

0.97 

0.54 

0.55 

0.54 

0.54 

SEVEM-Commander ... 

0.81 

0.83 

0.69 

0.67 

0.70 

0.70 

NILC-Commander .... 

0.85 

0.86 

0.64 

0.63 

0.66 

0.66 


is also performing very well in TTT and T+E, while the EEE tably, the fact that the clear degradation in the accuracy of the 
result is more than 2 (T/^l measurement for some methods does not seem to propa- 

The test described here has several limitations, the main and S^^e to T+E. 
most obvious one being that it has been performed on just two 

maps (due to lack of availability of a large sample of this type 7 4 Dependence on multipole number 
of simulation at this stage). Another clear issue is that some 

methods, in particular Commander, seem to perform much bet- In this section we discuss another stability test of our results, 
ter on actual data than on these simulations. On the other hand namely the dependence of the results for /nl on the maximum 
some important trends, observed in the data in previous tests, and minimum multipole number used in the analysis. This test 
are clearly reproduced here, like the good stability of SMICA and is most easily performed with the binned bispectrum estimator, 
SEVEN, especially for the combined T-\-E results and, most no- since it gets the dependence of /nl on £ for free with its stan- 
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Fig. 15. Scatter plots showing correlations between bispectrum 
modes extracted from the different Planck foreground-cleaned 
maps, for all possible pairs of component separation methods. 
Top: TTT bispectrum modes. Bottom: EEE bispectrum modes. 
While temperature shows a strong correlation, the loss of corre¬ 
lation in polarization between the different methods is evident in 
these plots, as discussed in the text and quantified in Tables 14 
and 15. Results here and in Table 15 have been obtained using 
the high resolution modal pipeline (2001 modes), while results 
in Table 14 have been obtained with the low resolution modal 
pipeline. By construction, the high resolution pipeline is measur¬ 
ing not the full EEE bispectrum of the map, but the component 
of EEE that is orthogonal to TTT. For this reason, measured by 
the two pipelines for EEE will not be identical. With this caveat 
in mind, the agreement between the two modal approaches is 
very good. 


dard analysis, simply by leaving out bins in the final sum when 
computing /nl (the binned equivalent of Eq. (39)). 

The dependence on ^max of the results for the three standard 
primordial shapes (local, equilateral, and orthogonal modes), is 
shown in Fig. 16, for T -only, £'-only, and full T+E. As mentioned 
in Sect. 6, the KSW and binned estimators use ^max = 2500 
for temperature, while the modal estimators use ^max = 2000. 
As can be seen in the figure, both the T-only and T+E results 
are basically unchanged between i = 2000 and i = 2500 for 


Table 15. The statistic, Eq.(58), showing the degree of corre¬ 
lation between measured bispectrum (3 coefficients for the com¬ 
ponent separation methods shown in Fig. 15 (upper three rows 
for temperature, lower for polarization). Correlation between 
SMICA, NILC, and SEVEN is excellent in temperature; however, it 
declines markedly for the latter in polarization. Note that results 
are from the high resolution Modal 2 pipeline using the orthog¬ 
onal EEE component only. 



SEVEM 

NILC 

Commander 

SMICA (T) 

0.95 

0.94 

0.63 

SEVEM (T) 


0.92 

0.66 

NILC (T) 



0.72 

SMICA (E) 

0.39 

0.89 

0.55 

SEVEM (E) 


0.30 

0.50 

NILC (E) 



0.43 


Table 16. Correlation coefficients between pairs of EEE bispec¬ 
trum modes, extracted using two different masks for each of the 
four component-separated maps. We compare correlations mea¬ 
sured from data with Monte Carlo averages over 100 Gaussian 
realizations. The simulations were processed through the differ¬ 
ent component separation pipelines in the same way as the data, 
but do not include any foreground component. According to this 
test, modal expansions are stable for a change of sky coverage, 
with measured correlations in full agreement with expectations 
from simulations. 


Method 

EEE 

Data Simul. 

SMICA . 

0.87 

0.87 

SEVEM . 

0.87 

0.87 

NILC. 

0.87 

0.87 

Commander . . . 

0.88 

0.87 


all three shapes, showing that this difference has no impact on 
the results (as was to be expected from the excellent agreement 
between estimators in Table 10). In fact, the values of /nl for T 
and T+E are reasonably stable (given their error bars) down to 
much lower values of ^max, of about 1000. 

On the polarization side we can draw a similar conclusion. 
The binned estimator uses ^max = 2000 for polarization, while 
the other estimators use ^max = 1500, but we see that results 
for E remain basically unchanged between i = 1500 and i = 
2000. Central values and error bars for E for all three shapes 
have clearly converged by ^ = 1500, and are in fact reasonably 
stable down to much lower values of about 700. 

As we noted in the 2013 analysis 
(Planck Collaboration XXIV 2014), when going to the much 
lower WMAP resolution of ^max - 500, we agree with the 
slightly high value of that the WMAP team reported 

(Bennett et al. 2013). This value is also confirmed when includ¬ 
ing polarization. One can clearly see the value of the higher 
resolution of Planck. 

The dependence on ^min is shown in Fig. 17. Here all esti¬ 
mators used the same values, . =2 and = 40. As ex- 
plained in Planck Collaboration VII (2016), not all systematic 
and foreground uncertainties in the low-^ HFI polarization data 
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Table 17. Comparison of component separation methods, using Gaussian FFP8 simulations. We firstly consider Gaussian, 
foreground-free simulations, with simulated noise for each frequency band, process them through each of the four foreground 
cleaning pipelines, and measure /nl for the three standard shapes (columns labelled with “Input map”). We then include fore¬ 
grounds and repeat the measurement, before applying the cleaning, and including realistic noise levels for each method (columns 
labelled with “Input map -r foregrounds”); this step is performed in order to get an idea of the level of contamination introduced by 
foregrounds, before cleaning. Finally, we apply the different component separation methods, and again estimate /nl from the final 
maps (columns labelled with “Cleaned map”). The discrepancies between /nl measured on the input map, and /nl extracted from 
the cleaned map, provide a figure of merit to assess how well foregrounds are subtracted by different methods. Results below have 
been obtained with the KSW estimator and the “cleaned map” results were also checked with the binned estimator. 



Input map 

Input map + foregrounds 

Cleaned map 

Local 

Equilateral Orthogonal 

Local Equilateral Orthogonal 

Local Equilateral Orthogonal 


SMICA 


T . 

5.2 + 

5.8 

29 + 

71 

-8 + 34 

-107.0 + 

5.8 

-23 + 

71 

27 + 34 

7.8 + 

5.8 

38 + 

71 

-20 + 34 

E . 

-39 + 

30 

-99 + 

133 

59 + 69 

-10 + 

30 

-154 + 

133 

-41 + 69 

-56 + 

30 

-120 + 

133 

65 + 34 

T+E . . . 

5.9 + 

5.1 

14 + 

45 

-20 + 22 

-118.0 + 

5.1 

-32 + 

45 

8 + 22 

8.3 + 

5.2 

14 + 

45 

-22 + 22 

SEVEM 
















T . 

5.6 + 

5.7 

32 + 

69 

-8 + 32 

-113.2 + 

5.7 

-8 + 

69 

34 + 32 

12.7 + 

5.7 

35 + 

69 

-25 + 32 

E . 

-17 + 

41 

-149 + 

175 

28 + 95 

-14 + 

41 

-171 + 

175 

-44 + 95 

-22 + 

41 

-120 + 

175 

41+95 

T+E . . . 

1.1 ± 

5.3 

12 + 

49 

-37 + 24 

-126.0 + 

5.3 

-29 + 

49 

-57 + 25 

13.0 + 

5.3 

11 + 

49 

-41 + 24 

NILC 
















T . 

5.1 + 

5.7 

32 + 

69 

-5 + 31 

-102.0 + 

5.7 

-14 + 

69 

32 + 31 

17.8 + 

5.7 

85 + 

69 

-16 + 31 

E . 

-52 + 

33 

-157 + 

156 

72 + 73 

-6 + 

33 

-155 + 

156 

-47 + 73 

-76 + 

33 

-179 + 

156 

113 + 73 

T+E . . . 

5.7 + 

5.0 

7 + 

46 

-15 + 21 

-117.0 + 

5.9 

-27 + 

46 

12 + 21 

15.8 + 

5.0 

-20 + 

46 

-7 + 21 

Commander 

T . 

0.5 + 

6.2 

-5 + 

73 

-14 + 36 

-127.0 + 

6.2 

-25 + 

73 

-137 + 36 

25.6 + 

6.2 

67 + 

73 

-17 + 36 

E . 

-51 + 

38 

-64 + 

160 

93 + 86 

-10 + 

38 

-153 + 

160 

-45 + 86 

-70 + 

38 

-78 + 

159 

138 + 86 

T+E . . . 

1.6 + 

5.4 

-2 + 

48 

-21 + 23 

-137.0 + 

5.4 

-29 + 

48 

13 + 23 

20.4 + 

5.4 

28 + 

48 

-11 + 23 


Table 18. Same test as in Table 17, but with a NG map as input, 
with = 8.8. For this case, we only report the final value 
after foreground cleaning for each method. Results below have 
been obtained with the KSW estimator and double-checked with 
the binned estimator. ISW-lensing contributions are removed. 


Cleaned map. Input /^/^^ = 8.8 


Local Equilateral Orthogonal 


SMICA 

T . 

3.1 

+ 

5.8 

47 + 

71 

-6 + 34 

E . 

-53 

+ 

30 

-113 + 

133 

94 + 69 

T+E . . . 

5.7 

+ 

5.1 

22 + 

45 

-19 + 22 

SEVEM 

T . 

8.0 

+ 

5.7 

43 + 

69 

-11 + 32 

E . 

-19 

+ 

41 

-112 + 

175 

35 + 95 

T+E . . . 

10.2 

+ 

5.3 

19 + 

49 

-37 + 24 

NILC 

T . 

10.2 

+ 

5.7 

84 + 

69 

7 + 31 

E . 

-76 

+ 

33 

-179 + 

156 

113 + 73 

T+E . . . 

10.1 

+ 

5.0 

20 + 

46 

4 + 21 

Commander 

T . 

22.2 

+ 

6.2 

81 + 

73 

-5 + 36 

E . 

-68 

+ 

38 

-78 + 

160 

132 + 86 

T+E . . . 

18.3 

+ 

5.4 

35 + 

48 

-9 + 23 


have been fully characterized yet, and hence it was decided to 
filter out these data. 

For equilateral and orthogonal shapes the values for /nl and 
their error bars are quite stable as a function of ^min up to about 


^ = 100 (and i ^ 300 for £'-only), which is not surprising, since 
these templates have little weight at low i. The local template, on 
the other hand, depends very strongly on the lowest multipoles, 
which is reflected in the very rapid growth of the error bars when 
^min increases. Looking at T-only and T+E we see a very similar 
pattern, with being reasonably stable, although there are 
some jumps. The local result for £'-only wanders a bit more out¬ 
side of the +1 (T region, in agreement with the other tests in this 
section, which also indicate that £’-only is not as stable as T -only 
and T+E. However, that is still quite acceptable, given the small 
weight of E'-only in the full T+E result. 

One can work out quite generally that when F is a subset 
of a data set X, and Px and Py are the values of a parameter 
P determined from these two data sets, then the variance of the 
difference Py - Px is equal to |Var(F 7 ) - Var(Fx)l Hence we can 
determine how likely the jumps in /nl as a function of ^min are. 
It turns out that the jump in the T-only value of between 
^min = 40 and^min = 53isa2.5(T effect (using the values of 
before subtraction of the ISW-lensing bias, which also depends 
on £). Similarly, the jump in the T-only value of between 
^min = 154 and ^min = 211 is (by chance) also a 2.5cr effect. 
Given the fact that there are 57 bins, having such a jump appears 
to be consistent from a statistical point of view. 


7.5. A directional analysis with a needlet-based modal 
estimator 

The validation tests on simulations, described in Sect. 7.3, point 
to SMICA and SEVEM as the best foreground cleaning methods 
for /nl estimation. Results in Table 10 show that SMICA also has 
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Fig. 16. Evolution of the /nl parameters (solid blue line with data points) and their uncertainties (dashed lines) for the three primor¬ 
dial bispectrum templates as a function of the maximum multipole number used in the analysis. From left to right the figures 
show, respectively local, equilateral, and orthogonal, while the different rows from top to bottom show results for T -only, £'-only, 
and full T+E. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of /nl (solid 
green lines without data points, showing the +1 cr range around the dashed green lines). The results are for SMICA, and assume all 
shapes to be independent, and that the ISW-lensing bias has been subtracted. They have been determined with the binned bispectrum 
estimator. 







slightly smaller error bars, thus making it the method of choice 
for our final results. 

As a further check of residual foreground contamination 
in the SMICA map, in this section we investigate the possi¬ 
ble directional dependence of SMICA-derived third-order statis¬ 
tics by means of a needlet-based modal estimator (i.e., an es¬ 
timator based on the decomposition described in Sect. 3.2, 
and references therein, where we use cubic combinations of 
needlets as our basis modes). In other words, we analyse the be¬ 
haviour of the needlet bispectrum (see Lan & Marinucci 2008, 
Rudjord et al. 2010 and Donzelli et al. 2012) on separate patches 
of the sky, and we study the fluctuations of the corresponding 
residuals. 

Rather than assuming a specific anisotropic model, we in¬ 
stead calculate the contribution to the local /nl from different 
regions of the sky and look for evidence of anisotropy in the 
result. 


Even though we focus here on directional contributions to the local 
shape (which was typically found as one of the most sensitive to residual 


Our modal needlet estimator has been validated with respect 
to the procedures considered in Sect. 5, showing excellent agree¬ 
ment. Since in this paper we use the needlet estimator only in this 
section, and as a diagnostic tool, we will not explicitly report the 
outcome of these validation tests here, for the sake of concise¬ 
ness. The advantages of using needlet-based modal estimators 
have been advocated in Lan & Marinucci (2008), Rudjord et al. 
(2010), Fergusson et al. (2010a), and Fergusson et al. (2012); we 
refer to these papers for more discussion and details. In short, 
however, they can be summarized as follows: 

1. it is possible to achieve a strong data compression, i.e., to 
investigate cubic statistics by means of a small number of 
modes (needlet frequencies) for many different bispectrum 
templates; 


foreground contamination), this type of directional analysis can be done 
in a model-independent way, by looking separately at different needlet 
modes; we leave this for future work. 


36 



















































Planck Collaboration: Planck 2015 Results. Constraints on primordial NG 



Fig. 17. Evolution of the /nl parameters (solid blue line with data points) and their uncertainties (dashed lines) for the three primor¬ 
dial bispectrum templates as a function of the minimum multipole number used in the analysis. From left to right the figures 
show respectively, local, equilateral, and orthogonal, while the different rows from top to bottom show results for T -only, £'-only, 
and full T+E. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of /nl (solid 
green lines without data points, with +1 cr range around the dashed green line). The results are for SMICA, and assume all shapes to 
be independent, as well as that the ISW-lensing bias has been subtracted. They have been determined with the binned bispectrum 
estimator. 


2. needlet transforms have good correlation properties in pixel 
space, which allows study of the pixel contribution to the 
/nl signal for the templates under investigation, by treating 
different directions independently. 

In our analysis, we first divide the sky into several large “re¬ 
gions”, with boundaries defined by the pixels of a HEALPix grid 
(Gorski et al. 2005). at lower resolution than the starting map 
(which is at A^side = 2048). For the low resolution grid we con¬ 
sider A^side values between 2 and 8. For each pixel in the coarse 
grid, we then compute the local /nl using our modal needlet es¬ 
timator, and neglect contributions from external regions. 

The correlation matrices between /nl measurements in dif¬ 
ferent regions were computed via Monte Carlo simulations, 
^side = 4 was chosen, as providing the best tradeoff between 
having a large number of regions for directional analysis (i.e., 
the total number of pixels in the low resolution grid), and a low 
correlation between different regions. This is shown in Fig. 18. 
It is readily seen that, at A^side = 4, the correlation is largely con¬ 
centrated in one or two points near the main diagonals (where it 


is still low, never exceeding 34 %), and falls off rapidly for all 
other pixels. Note that the results here refer to temperature only. 
The EEE local polarization error bars are large, even for the full 
sky analysis, making this directional approach uninformative for 
Planck polarization data. We concentrate on the TTT bispectrum 
here, complementing other validation tests in this section, which 
are mostly focused on polarization. 

Having obtained our correlation matrices, and having shown 
that different regions are essentially uncorrelated, we can then 
proceed to extract /nl for each region in the actual SMICA map. 
As a test of directional-dependent contamination of the /nl mea¬ 
surement, we can also compare our results, region by region, 
with the fluctuations expected by looking at the diagonal of our 
Monte Carlo correlation matrix. 

The results of this analysis are shown in Fig. 19. In the left 
panel, we represent the directional local /nl map, extracted with 
this method. In the right panel, we report the /nl values, region 
by region, and compare them to expectations from simulations. 
The red line gives the expected standard deviation, while the blue 
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Pixel number (Healpix RING scheme) 


Fig. 18. Temperature only, local shape, pixel correlation matrix 
from Monte Carlo analysis, at A^side = 4. The most correlated 
pixels are those closest to the main diagonal; however, these val¬ 
ues are always lower than 34 % in the chosen case. 

one gives estimates on the component-separated maps with the 
Monte Carlo error bars. Our estimator is normalized in such a 
way that the sum of all these contributions would yield exactly 
the /nl estimator for the full map, so these results can be viewed 
as a partition of the estimates along the different directions. It is 
readily seen that no significant fluctuation occurs, so that our re¬ 
sults are consistent with the absence of directionally-dependent 
features (which could occur due to, for instance, residual fore¬ 
ground contamination). As an additional check, we have also 
investigated the possible presence of a dipole in these data, and 
found that our results are consistent with Gaussian isotropic sim¬ 
ulations. 

7.6. Summary of the main validation results 

Throughout Sect. 7 we have shown a battery of tests aimed at 
evaluating the robustness of our data set, from the point of view 
of bispectrum estimators, focusing especially on the polarizaton 
part. We studied the stability of our results (local, equilateral, 
and orthogonal /nl measurements, model-independent bispec¬ 
trum reconstruction) under a change of sky coverage, multipole 
range, and choice of component separation methods. We also 
considered simulated data sets and studied the ability of differ¬ 
ent component separation methods to recover the input /nl after 
foreground subtraction. Our main conclusions from these tests 
can be summarized as follows: 

“ TTT and T+E results are stable both in the pixel and har¬ 
monic domains, for different component separation meth¬ 
ods. For SMICA, we also checked that TTT temperature con¬ 
straints on local /nl show no evidence of a directional varia¬ 
tion via a needlet-based analysis. 

- SMICA and SEVEM perform better than NILC and Commander 
at recovering the original /nl in foreground-cleaned simula¬ 
tions. At the same time, SMICA allows slightly better con¬ 
straints on /nl than SEVEM, due to a lower (by a small 
amount) noise level. 

- EEE bispectra, and related /nl measurements, have some 
problems, and do not pass all the tests. Different compo¬ 
nent separation methods show a low level of consistency 
(especially when comparing pixel-based cleaning methods 
to harmonic-based cleaning methods). This disagreement is 


only partly alleviated by choosing a larger £’-mask, so that 
residual foregrounds do not seem to fully explain all issues. 
An important caveat, already pointed out previously, is that 
the noise model in polarized FFP8 simulations is known to 
underestimate the actual noise level in the data, leading to 
some degree of underestimation of the error bars in our EEE 
results. We stress again, however, that this has little impact 
on the final T+E constraints, due to the high noise level and 
consequent low statistical weight of EEE bispectra. This was 
verified in detail both on data and on simulations. 

In light of the above analysis, we conclude that, as far as 
bispectrum estimation is concerned, our best cleaned map is the 
one produced by SMICA, in line with our previous 2013 analysis. 
We also conclude that our rTT-based /nl constraints, summa¬ 
rized in Tables 10 and 11, are robust. Joint T+E constraints pass 
all our validation tests. On the other hand, in the light of the re¬ 
maining issues in the EEE bispectra and in the FFP8 polarized 
simulations, as we stressed at the end of Sect. 6.1, we suggest 
that all measurements that include polarization data in this paper 
should be regarded as preliminary. 

8. Other non-Gaussianity shapes for /nl 

This section discusses new searches for NG beyond standard 
single-field inflation. The focus here is on extensions to the 
analysis undertaken in Planck Collaboration XXIV (2014) with 
new limits on isocurvature models, further oscillatory models 
over a broader frequency range, and parity-violating tensor NG. 
However, we also briefly revisit all the non-standard models con¬ 
strained in Planck Collaboration XXIV (2014), including effec¬ 
tive field theory, non-Bunch Davies, and directionally-dependent 
models, in particular noting the impact of the new (preliminary) 
polarization data on the previous constraints. 

8.1. Isocurvature non-Gaussianity 

Here we show the results obtained for a study of the isocurva¬ 
ture NG in the Planck 2015 SMICA map using the binned bis¬ 
pectrum estimator. As explained in Sect. 2.3, we only investi¬ 
gate isocurvature NG of the local type, and in addition always 
consider only one isocurvature mode (either cold dark matter, 
neutrino density, or neutrino velocity isocurvature) in addition 
to the adiabatic mode. Note that the baryon isocurvature mode 
behaves identically to the cold dark matter one, only rescaled by 
factors of Db/^c. so there is no need to consider it separately. 
In that case there are six different /nl parameters: a purely adia¬ 
batic one (a,aa), which corresponds to the result from Sect. 6), a 
purely isocurvature one (i,ii); and four mixed ones (see Sect. 2.3 
for an explanation of the notation). 

The results are given in Table 19^^ and show no clear signs 
of any isocurvature NG. There are a few values that deviate from 
zero by up to about 2.5 cr, but such a small deviation, in partic¬ 
ular given the large number of results, cannot be considered a 
detection. We do see that many constraints are tightened consid¬ 
erably when including polarization, by up to the predicted factor 
of about six for the cold dark matter a,ii, i,ai, and i,ii modes in 

Compared to definitions in the literature based on ^ and S (see 
e.g., Langlois & van Tent 2012), here we adopt definitions based on 
^adi = 3^/5 and Oiso = 5'/5, in order to make the link with the standard 
adiabatic result more direct. Conversion factors to obtain results based 
on ^ and S are 6/5, 2/5, 2/15, 18/5, 6/5, and 2/5, for the six modes, 
respectively. 


38 













Planck Collaboration: Planck 2015 Results. Constraints on primordial NG 



Adimensional 



Fig. 19. Temperature only, local /nl directional contributions from SMICA. As explained in the text, summing over all the pixel 
values would give the full sky /nl needlet estimator result. The left panel displays the directional /nl map. On the right, the blue 
points represent the /nl contibution for each direction (i.e., for each pixel in the directional map), with Monte Carlo error bars. The 
red line is the average from simulations, which is consistent with zero. 


Table 19. Results for local isocurvature NG, determined from the SMICA Planck 2015 map with the binned bispectrum estimator. In 
each case the adiabatic mode is considered together with one isocurvature mode (either cold dark matter, neutrino density, or neutrino 
velocity isocurvature). As explained in the text this gives six different /nl parameters, indicated by the different combinations of 
the adiabatic (a) and isocurvature (i) modes. Results with two significant digits are shown for both an independent and a fully joint 
analysis, for T-only, £’-only, and full T+E. In all cases the ISW-lensing bias has been subtracted. 


/nl 


Independent Joint 

Shape Cold dark matter Neutrino density Neutrino velocity Cold dark matter Neutrino density Neutrino velocity 


T a,aa. 

1.3 

+ 

5.4 

1.3 

+ 

5.4 

1.3 

+ 

5.4 

21 + 

13 

-27 + 

52 

-32 + 

48 

T a,ai. 

-2 

+ 

10 

-4 

+ 

15 

47 

+ 

29 

-39 + 

26 

140 + 

210 

370 + 

350 

T a,ii. 

59 

+ 

910 

-130 

+ 

280 

750 

+ 

360 

17000 + 

8200 

-4500 + 

4500 

-1300 + 

3800 

T i,aa. 

6 

+ 

50 

3.0 

+ 

9.0 

1.0 

+ 

4.7 

96 + 

120 

40 + 

99 

-27 + 

51 

T i,ai. 

3 

+ 

66 

-5 

+ 

22 

26 

+ 

21 

-2100 + 

1000 

220 + 

630 

75 + 

170 

T i,ii . 

76 

+ 

280 

-100 

+ 

250 

440 

+ 

230 

4200 + 

2000 

-750 + 

2400 

-970 + 

1400 

E a,aa. 

34 

+ 

34 

34 

+ 

34 

34 

+ 

34 

66 + 

50 

51 + 

120 

-140 + 

150 

E a,ai. 

-31 

+ 

200 

70 

+ 

140 

78 

+ 

93 

-380 + 

310 

-280 + 

640 

1100 + 

620 

E a,ii. 

-4200 

+ 

4000 

-520 

+ 

2300 

190 

+ 

940 

-8800 + 

6100 

-6400 + 

6200 

-9400 + 

3900 

E i,aa. 

-10 

+ 

87 

42 

+ 

42 

23 

+ 

27 

27 + 

180 

52 + 

170 

54 + 

120 

E i,ai. 

94 

+ 

250 

83 

+ 

130 

45 

+ 

62 

910 + 

770 

670 + 

850 

-190 + 

420 

E i,ii . 

690 

+ 

2200 

390 

+ 

1400 

260 

+ 

460 

-6000 + 

5300 

-4100 + 

5300 

2200 + 

1600 

T+E a,aa .... 

0.7 

+ 

4.9 

0.7 

+ 

4.9 

0.7 

+ 

4.9 

5 + 

10 

-35 + 

27 

2 + 

24 

T+E a,ai .... 

-2.6 

+ 

9.7 

-5 

+ 

14 

17 

+ 

22 

-12 + 

20 

74 + 

94 

330 + 

130 

T+E a,ii. 

130 

+ 

450 

-130 

+ 

240 

130 

+ 

230 

-1800 + 

1300 

-3000 + 

1400 

-3200 + 

1200 

T+E i,aa .... 

30 

+ 

26 

5.6 

+ 

7.7 

-0.7 

+ 

4.1 

53 + 

47 

51 + 

45 

-44 + 

24 

T+E i,ai. 

26 

+ 

38 

2 

+ 

19 

6 

+ 

15 

140 + 

170 

170 + 

210 

20 + 

74 

T+E i,ii. 

38 

+ 

170 

-26 

+ 

180 

85 

+ 

130 

-280 + 

390 

-390 + 

860 

480 + 

430 


the joint analysis. As discussed in detail in Sect. 7, results in¬ 
cluding polarization data should be considered preliminary, and 
that is even more important here, since these results depend so 
strongly on the additional information from polarization. 

In the results so far we allowed for a possible correlation 
between the isocurvature and adiabatic modes. However, if we 
assume that they are completely uncorrelated, with a zero cross 
power spectrum, then there are only two /nl parameters, the a,aa 
and i,ii ones. In Table 20 we give the results for this uncorrelated 


case. The independent results are the same as in the previous 
table, while in the joint results one can clearly see the difference 
between the neutrino density mode (the bispectrum template of 
which has a large overlap with the adiabatic one), and the cold 
dark matter and neutrino velocity modes (with templates that are 
very different from the adiabatic one). Again there is no evidence 
for any isocurvature NG; the almost 2 cr result for the neutrino 
velocity isocurvature mode in the temperature-only case does 
not survive the addition of polarization. 
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Table 20. Similar to Table 19, except that we now assume that the adiabatic and isocurvature mode are completely uncorrelated. 
Hence there are only two /nl parameters in this case, a purely adiabatic one and a purely isocurvature one. 


/nl 

Independent Joint 

Shape Cold dark matter Neutrino density Neutrino velocity Cold dark matter Neutrino density Neutrino velocity 

ra,aa. 1.3+ 5.4 1.3+ 5.4 1.3+ 5.4 1.0+ 5.3 19 + 12 -0.2+ 5.4 

ri,ii. 76 + 280 -100 + 250 440 +230 65 + 280 -840 + 540 440 +230 

E a,aa. 34 + 34 34 + 34 34 + 34 33 + 35 42 + 40 35 + 40 

£i,ii. 690 +2200 390 + 1400 260 +460 210 +2200 -680 + 1700 -31 +540 

T+E a,aa_ 0.7 + 4.9 0.7 + 4.9 0.7 + 4.9 0.5 + 5.0 3.0 + 7.9 -0.3 + 4.9 

T+EiM . 38 + 170 -26 + 180 85 + 130 35 + 170 -120 + 290 87 + 130 


8.2. Feature models 

An important and well-motivated class of scale-dependent bis¬ 
pectra is the feature model, characterized by linear oscillations 
described by Eq. (15) and its variants in Eqs. (16) and (18). 
In Planck Collaboration XXIV (2014) we performed an initial 
search for a variety of feature models using the Modal estimator. 
This earlier search was limited to oj < 200 by the native reso¬ 
lution of our implementation of the Modal estimator (using 600 
modes), roughly the same range as the initial WMAP bispectrum 
feature model searches at lower precision (with only 50 eigen- 
modes (Eergusson et al. 2012). Note, in the previous Planck re¬ 
lease we used wavenumber k^ in line with the theory literature, 
but here we switch to frequency oj, in line with more recent 
observational power spectrum searches; the two are related by 
OJ = 27113kc. With the improved estimator resolution (now using 
2001 modes) we are able to achieve convergence over a broader 
range up to tu = 350. We perform a frequency scan of 350 sam¬ 
pling points between oj = 10 and oj = 350, i.e., 35 independent 
frequencies and 10 phases. We also extend the number and vari¬ 
ety of feature and resonance models that are investigated, essen¬ 
tially probing the resolution domain in which we have obtained 
a reliable Modal bispectrum reconstruction (see Eig. 4). 

Constant feature ansatz: Eor the constant feature shape of 
Eq. (15), we can extend the frequency range much further with 
another approach. As the bispectrum in Eq. (15) is separable, it 
allows the construction of a KSW estimator (Miinchmeyer et al. 
2014) for direct bispectrum estimation at any given frequency. 
The bispectrum can be written as a sum of sine and cosine com¬ 
ponents which can be estimated separately (equivalent to mea¬ 
suring the amplitude and phase above) and this method was used 
to constrain frequencies up to tu = 3000. The range where the 
two estimators overlap provides validation of the two methods 
and excellent agreement was seen (see Eig. 20). 

Apart from cross-validation with two estimators, we have 
undertaken further tests to determine the robustness of the re¬ 
sults to foreground and noise effects. In Eig. 21 (left panel), we 
show the effect on feature model results of the subtraction of 
the simple point source bispectrum, as well as the ISW-lensing 
bispectrum. This study was a major motivation for adopting the 
more conservative “extended” common mask, because the con¬ 
sistency between different component separation methods im¬ 
proved markedly for low frequencies, with the original common 
mask requiring much larger point-source subtractions (e.g., for 
NILC subtraction the maximum raw significance has reduced 
from cr = 4.0 to (Xciean = 2.2 dit OJ = 110). After cleaning 
these signals, the SMICA, SEVEM and NILC results are in good 


agreement, and also consistent with each other when polariza¬ 
tion is included (while the Commander results generally have 
a larger variance and so are not included in the plotted av¬ 
erages). Eortunately, the effect of subtracting ISW-lensing and 
point source bispectra diminishes rapidly at higher frequencies 
OJ > 200 and should be negligible; subtraction was only under¬ 
taken in the Modal region i < 350. In Eig. 21 (right panel), 
we show the effect on the averaged significance of reducing the 
Planck domain from the usual T’^ax = 2000 to T’max = 1500 
(^max = 1500 to T’max = 1000 for £’-modes). Despite the non¬ 
trivial change in overall signal-to-noise entailed, there is no 
strong evidence for an T’-dependent signal, as might be expected 
if there was substantial NG feature contamination in the noise- 
dominated region. Einally, we note that most peaks at low oj 
show some correlation between T-only and T+E, although there 
are notable exceptions, such as the peak dXoj = 180, which is re¬ 
moved after inclusion of polarization (see also the phase plots in 
Eig. 20 before marginalization). The temperature feature peaks 
observed in Eig. 21 at tu 110, 150, and 180 are consistent with 
the peaks identified previously in Planck Collaboration XXIV 
(2014). 

In Eig. 20, the full set of frequency results, 0 < tu < 3000, 
for the constant feature model in Eq. (15) is shown for both the 
Modal and KSW estimators. There is good agreement for re¬ 
gions where the estimators overlap and PS and ISW lensing bis¬ 
pectra subtraction is not necessary {oj > 200). Generally, there 
is tighter consistency between temperature-only results, than 
with polarization where there is additional scatter between fore¬ 
ground separation methods. Scanning across the full frequency 
range, there is no strong evidence for any large maximum that 
might indicate unequivocal evidence for a feature model signal. 
The maximum peak significance obtained with either T -only or 
T+E is consistent with expectations for a Gaussian model over 
this frequency range. In particular, for the KSW estimator us¬ 
ing the SMICA data, the highest significance found in the range 
200 < OJ < 3000 is 3.2 cr in T-only and 2.9 cr in T+E. To 
gauge the likelihood of these results occurring randomly, real¬ 
istic Gaussian SMICA simulations were analysed with the KSW 
estimator and found to typically produce a highest peak with 
3.1 (+0.3) cr over the same frequency range. 

Generalized feature models: We have also deployed the Modal 
estimator to look at (non-separable) feature models with equilat¬ 
eral and flattened cross-sections, as motivated by varying sound 
speed scenarios and those with highly excited states, respec¬ 
tively. In the left panels in Pig. 22 we show results from the equi¬ 
lateral feature model of Eq. (16), including the frequency/phase 
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Fig. 20. Constant feature model results for both T-only and T+E data across a wide frequency range. The upper four panels show the feature 
signal in the Modal range 0 < du < 350. The two upper left panels show contours of the raw significance cr obtained from the SMICA map as a 
function of the frequency du, for T -only and T+E, respectively. The upper right panels show the maximum signal after marginalizing over phase 0 
for both T-only and T+E for all foreground separation models. The third and fourth panels show the maximum feature signal in both T-only and 
T+E across the frequency range 0 < du < 1000, plotting both Modal results (dashed lines) and KSW results (solid lines for 200 < oj < 1000); these 
show good agreement in the overlap. The lower two panels give the maximum KSW results for T-only and T+E in the range 1000 < oj < 3000). 
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Fig. 21. Constant feature ansatz validation for the Modal estimator, showing the effect of ISW-lensing and point source subtraction for i < 300 
(left panel) and the impact of a lower I maximum cutoff on the average signal (right panel), i.e., lowering ^max = 2000 to 1500 {T) and ^max = 1500 
to 1000 {E). All Modal 2 results in Sect. 8 have used the extended common mask, except the validation analysis at different resolutions (right 
panel) which for consistency employs the common mask. 



Fig. 22. Generalized feature models analysed at ^max = 2000 (F-modes ^max = 1500) for the different Planck foreground separation methods, 
SMICA (blue), SEVEM (red), NILC (green), and Commander (yellow), together with the SSN (SMICA - SEVEM - NILC ) average (black). The left three 
panels apply to the equilateral feature models, showing, respectively, in the top panel the full feature survey significance at each frequency and 
phase (temperature only), the maximum significance at each frequency for temperature only (middle), and with polarization (lower). The right 
three panels apply to the flattened feature models, showing the maximum significance at each frequency for temperature only (top right) and with 
polarization (middle right), along with significance at each frequency and phase for temperature and polarization (right lower). 


contours before marginalization for the SMICA T-only map. 
Multiple peaks are apparent in the temperature signal across 
the Modal range up to an average maximum 3.3 cr raw signif¬ 
icance. However, from the lower panel it is clear that the po¬ 
larization signal is not well correlated with the temperature in 
the equilateral case, generally reducing peak heights with the 
maximum now about 2.6 cr (while eliminating the tu = 180 
peak altogether). This temperature peak remains present with 
the ^niax = 1500 cutoff, where the signal is slightly higher, 
but the polarization in this case is less well correlated (using 
^max - 1000)- For Gaussian noise we would not expect polar¬ 
ization to reinforce a high temperature signal on average. It may 
also be that the equilateral temperature signal has some resid¬ 


ual diffuse point-source contamination. The equilateral feature 
model is the most affected by the removal of point sources, so 
the presence of a more complex correlated PS bispectrum (not 
removed by the constant PS template subtraction) remains for 
future investigation. 

Results for the flattened feature model in Eq. (18) are shown 
in Fig. 22 (right panels), displaying more coherence between 
temperature and polarization. The temperature signal with a 
2.6 cr peak between 50 < tu < 150 is reinforced by polariza¬ 
tion and merges to make a broad 3 a peak around oj = 90, to¬ 
gether with another distinct peak Sit oj ^ 140. Such broad fre¬ 
quency peaks are not expected, because neighbouring feature 
models should be nearly uncorrelated over a range AtUeft ~ 13 
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Fig. 23. Single field feature model significance with a cos cjjK scaling dependence (Eq. 19) (left panels, T-only upper and T+E lower) or 
with a K sin coK scaling (Eq. 20) (right panels). To find the maximum signal, these results have been marginalized over the cr-dependent envelope 
function ranging from cr ^ 0 (no envelope) to the maximum cutoff allowed by the Modal resolution acj - 90. 

Table 21. Peak statistics for the different feature models showing the Raw peak maximum significance (for the given Modal sur¬ 
vey domain), the corrected significance of this Single maximum peak after accounting for the parameter survey size (the “look- 
elsewhere” effect) and the Multi-pedk statistic which integrates across the adjusted significance of all peaks to determine con¬ 
sistency with Gaussianity. SMICA, SEVEN and NILC map analyses exhibit satisfactory bispectrum agreement for all the different 
models, whereas the Commander results produce some anomalously large results, especially for polarization. The significant signal 
for the equilateral features model in the T -only multi-peak statistic is reduced when polarization is added. The flattened feature 
model produces interesting results, which are reinforced with polarization to the 3 cr-level, with a high multi-peak significance. 




SMICA 



SEVEM 



NILC 


Commander 


Raw 

Single 

Multi 

Raw 

Single 

Multi 

Raw 

Single 

Multi 

Raw 

Single 

Multi 

Features constant T-only 

2.7 

0.5 

1.7 

2.6 

0.4 

1.5 

2.8 

0.7 

2.2 

2.9 

0.8 

2.7 

Features constant T+E 

2.7 

0.5 

1.9 

2.8 

0.7 

2.5 

2.8 

0.7 

2.4 

2.6 

0.4 

1.5 

Features equilateral T -only 

3.3 

1.5 

4.0 

3.2 

1.3 

3.5 

3.3 

1.6 

4.1 

2.9 

0.9 

2.5 

Features equilateral T+E 

2.6 

0.4 

1.3 

2.6 

0.4 

1.6 

2.8 

0.7 

1.9 

2.7 

0.6 

1.5 

Features flattened T -only 

2.5 

0.3 

1.4 

2.6 

0.4 

1.6 

2.7 

0.5 

2.1 

2.8 

0.8 

2.7 

Features flattened T+E 

2.9 

0.9 

2.9 

3.0 

1.1 

3.5 

3.1 

1.2 

3.8 

3.1 

1.2 

3.8 

cos features T -only 

2.5 

0.7 

1.9 

2.3 

0.6 

1.6 

2.7 

1.0 

2.5 

2.2 

0.3 

1.1 

cos features T+E 

2.7 

1.0 

2.5 

2.7 

1.1 

2.6 

2.6 

1.0 

2.5 

2.4 

0.6 

1.8 

K sin feature T -only 

2.8 

1.2 

2.8 

2.7 

1.1 

2.7 

3.0 

1.5 

3.4 

2.6 

0.9 

2.3 

K sin features T+E 

2.1 

0.3 

1.0 

2.9 

1.4 

3.1 

2.4 

0.6 

1.7 

2.3 

0.5 

1.6 


(as we discuss below). As the phase plots in Fig. (22) indicate, 
this breadth in frequency oj may reflect the neighbouring feature 
models adjusting phase 0 to match an underlying NG signal of 
a related, but different, nature. We note also that the frequency 
region for oj < 100 is susceptible to some degeneracy with cos¬ 
mological parameters. We shall consider a “look-elsewhere” sta¬ 
tistical analysis of these results below. 

Single-field feature solutions: We have also searched for the spe¬ 
cific analytic solutions predicted for single-field inflation models 
with step-like potential features, as given in Eqs. (19) and (20), 
with results shown in Fig. 23. The highest peaks for the K^- 
cosine model occur around 2.5 a with temperature-only, then 
rises to 2.7 cr when polarization data are included, again with 
peaks at other distinct frequencies apparent. The K-sino model 
shows a similar apparent signal level, with a maximum T-only 
2.7 (T peak, dropping to 2.4 cr with polarization. One further 
difficulty with a positive interpretation of these bispectrum re¬ 
sults in this low frequency range is that stronger S/N counter¬ 
parts in the power spectrum are predicted for these simple mod¬ 


els (Adshead et al. 2011), whereas no significant correlated os¬ 
cillation signals are apparent at the relevant peak frequencies 
Planck Collaboration XX (2016). 

Feature model peak statistics: In order to determine consis¬ 
tency with Gaussianity for these feature model results, we can 
apply a number of statistical tests developed for this specific 
purpose (Fergusson et al. 2015) and, if warranted, also apply 
these jointly in combination with power spectrum results, as for 
the WMAP polyspectra analysis (Fergusson et al. 2014). When 
scanning across the (tu, 0) parameter-dependent feature mod¬ 
els, we are searching through independent models for which 
Gaussian noise, by chance, can lead to a large apparent signal. 
We must correct for this multiplicity of tests when determining 
the actual significance of results for a given model — this is a 
quantitative correction for any model with free parameters, dis¬ 
tinct from the a posteriori choice of models to test. The simplest 
approach is to determine whether the maximum peak is consis¬ 
tent with Gaussian expectations, which can be determined from 
Monte Carlo simulations. However, in Fergusson et al. (2015) 
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it was recognized that these feature models can be accurately 
characterized analytically with a ;^-distribution having two de¬ 
grees of freedom^^. Taking = 2000 and using this analysis, 
the frequency step size between models that are uncorrelated is 
approximately AtUeff = 13, so we have an effective number of 
independent feature models A^eff ~ 27 for the Modal frequency 
range (with A^eff ~ 230 across the larger constant feature sur¬ 
vey range). Accordingly appropriate look-elsewhere corrections 
have been applied to find an adjusted significance for the max¬ 
imum peak signal found in all the feature model searches un¬ 
dertaken, which is shown in Table 21. Given that this feature 
model survey is over many independent frequency models and 
combinations of data, even the highest raw significances above 
3 cr (“Raw” column in Table 21) are reduced to a corrected sig¬ 
nificance below 2(T (“Single” column). Hence, there appears to 
be no evidence from maximum peak statistics for feature model 
deviations from Gaussianity. 

Nevertheless, we also examine the possibility that multiple 
feature models are contributing to a NG signal, given the appar¬ 
ent emergence of several preferred frequency peaks. This inte¬ 
grated multi-peak statistic can also be accurately approximated 
analytically (Eq. 60) using a ;^-distribution; essentially we sum 
over all independent frequencies using the single peak signif¬ 
icance adjusted for the “look-elsewhere” effect (see Eq. (59)). 
Most of the signal surveys exhibit an unusual number of broad 
overlapping peaks within the accessible frequency domain, so 
the multi-peak statistic does yield a much higher significance, 
with many models above 2 cr after “look-elsewhere” correction. 
Notable cases are the temperature-only signal for the equilat¬ 
eral feature model which yields an average significance of 3.4 cr 
across the foreground-cleaned maps with concordant bispectrum 
results (i.e., SMICA, SEVEN and NILC); however, this interesting 
multi-peak significance drops to only 1.6 cr when the polariza¬ 
tion data are included (assuming the reliability of E results). On 
the other hand, the flattened feature model has an average multi¬ 
peak significance of 1.7 cr in temperature only, which rises to 
3.4 cr with polarization included (higher at 3.7 cr if Coimnander 
data were to be included in the average). In this case, beyond 
the number of peaks, it is also their width that contributes, with 
the main signal around oj ^ 90, much broader than Atu 13. 
Einally, after look-elsewhere effects are taken into account, the 
-cosine single-field solutions yield multi-peak statistics that 
rise with the inclusion of polarization data from 2.0 cr to 2.5 cr, 
while the Af-sine falls from 3.0 cr (T) to 1.9 cr {T-\-E). 


For the feature model with parameters (o, 0), the adjusted sig¬ 
nificance S for the raw significance cr after accounting for the “look- 
elsewhere” effect is given by (Fergusson et al. 2015) 

5 = V2Err‘[(F,,,2(o-))'^'f], (59) 

where is the cumulative distribution function of the ;^f-distribution 
of degree two and Aeff is the effective number of independent feature 
models. This can also be used to investigate whether feature models are 
contributing at several frequencies. This multi-peak statistic integrates 
over all peak signals using the corrected significance 5, i.e., 

S"=^y] , (60) 

AWeff ^ 

where Adu is the sampling step-size and AdUeff is the effective corre¬ 
lation scale between independent models given by AdUeff = (Omax - 
c>>min)/(A^eff- !)• Esscntially this sums up all significant peaks that have a 
non-zero adjusted significance, after accounting for the look-elsewhere 
effect in Eq. (59). 


An interesting, but not entirely coherent, picture emerges 
from these searches for non-standard models in the new 
Planck temperature data, especially when combined with 
the additional (preliminary) polarization information. In 
Planck Collaboration XXIV (2014), we noted that the feature 
model searches provided interesting hints of NG. This more rig¬ 
orous statistical analysis confirms this view, allowing for sev¬ 
eral alternative feature model explanations of the apparently 
high NG signal observed in the bispectrum reconstructions (see 
Eig. 4). However, there is no strong evidence for a single large 
feature model at a particular frequency; rather, the high multi¬ 
peak statistics indicate signal that is spread across several broad 
peaks. Given the variability between different feature models 
and polarization component-separation methods, we note the 
caveat that the integrated multi-peak statistic could be sensitive 
to calibration issues and foreground contamination. Eor this rea¬ 
son, we do not make strong claims for these non-standard signals 
at this stage, but we note that oscillatory models will continue to 
be investigated thoroughly over a wider frequency domain and 
using the more reliable polarization data available in the final 
Planck data release. 

8.3. Resonance/axion monodromy nnodels 

Generalized resonance models: Using the Modal expansion, we 
have embarked on a survey of the simplest resonance model 
(Eq. 10), as well as the equilateral and flattened variants pro¬ 
posed in the literature, i.e., described by Eqs. (13) and (14), re¬ 
spectively. The raw significance for the resonance models for 
both temperature-only and temperature plus polarization data 
are shown in Eig. 24; these are the maximal results marginalized 
over the phase parameter 0. Previously, the resonance model was 
studied in Planck Collaboration XXIV (2014) using the Modal 
expansion over a narrower frequency range yielding no results 
above a raw significance of 1 cr. In this extended analysis over 
a wider frequency range, the constant sin(log) model (Eq. 10) 
produces 2.2 cr peaks for T-only, and 2.6 cr for T+E. The equi¬ 
lateral resonance model (Eq. 13) achieves a maximum 2.8cr in 
T-only diioj ^ 35, rising to a more impressive average 3.2 cr for 
TfE. Eor the flattened case (Eq. 14) we have 2.5 cr and 3.0 cr, re¬ 
spectively diioj ^ 12. Qualitatively, the results shown in Eig. 24, 
exhibiting broad peaks, are similar to those for feature models. 

Resonance model peak statistics: To determine the statistical sig¬ 
nificance of these results given the look-elsewhere effect of scan¬ 
ning across the parameters (tu, 0), we have used the two peak 
statistics defined above in Eqs. (59) and (60) for feature models 
(Eergusson et al. 2015). In this case, the maximum peak statis¬ 
tic for the constant resonance model of 2.6 cr (T+E) is readjusted 
to an unremarkable “look-elsewhere” single peak significance of 
0.9 cr. Likewise the apparently significant results above 3 cr for 
the equilateral and flattened models now fall below 2.0 cr with 
TfE. Using the single peak statistic alone, we would conclude 
that there is no strong evidence for any individual resonance 
model. Resonance models also generate oscillations in the power 
spectrum, and an analysis based on the 2015 temperature and po¬ 
larization likelihood is presented in (Planck Collaboration XX 
2016). A combined statistical treatment of resonance model 
power spectrum and bispectrum results will be reported in the 
future. 

The multi-peak statistic in Eq. (60) integrates the resonance 
model signal across all frequencies to determine consistency 
with Gaussianity. The constant resonance model has a modest 
multi-peak signal but, like the feature models, the equilateral and 
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Fig. 24. Generalized resonance models analysed at ^max = 2000 (F-modes ^^ax = 1500) for the different Planck foreground separation methods, 
SMICA (blue), SEVEM (red), NILC (green), and Commander (yellow), together with the SSN (SMICA - SEVEM - NILC ) average (black). The 
upper panels apply to the constant resonance model (Eq. 10), with T-only (left) and T+E (right), the middle panels give results for the equilateral 
resonance model (Eq. 13), and the lower panels for the flattened resonance model (Eq. 14). Both the equilateral and flattened resonance models 
produce broad peaks, which are reinforced with polarization (middle and bottom right panels). 

Table 22. Peak statistics for the resonance models showing the maximum Raw peak significance, the Single peak significance 
after accounting for the parameter survey “look-elsewhere” effect, and the Multi-peak statistic integrating across all peaks (also 
accounting for the “look-elsewhere” correction). There is some evidence for a high signal for both the equilateral and flattened 
resonance models, which increases when the polarization signal is added. This table does not include the results of the high frequency 
resonance model estimator, whose significance was assessed independently. 




SMICA 



SEVEM 



NILC 



Commander 



Raw 

Single 

Multi 

Raw 

Single 

Multi 

Raw 

Single 

Multi 

Max. 

Single 

Multi 

Sin(log) constant T -only 

2.4 

0.7 

1.2 

2.2 

0.4 

0.9 

2.0 

0.2 

0.7 

2.5 

0.8 

1.6 

Sin(log) constant T+E 

2.4 

0.7 

1.7 

2.7 

1.1 

2.4 

2.6 

1.0 

2.2 

2.7 

1.1 

2.5 

Sin(log) equilateral T -only 

3.0 

1.6 

2.4 

2.8 

1.2 

2.0 

2.5 

0.9 

1.5 

2.6 

1.0 

2.1 

Sin(log) equilateral T+E 

3.5 

2.2 

3.5 

3.0 

1.5 

2.8 

3.1 

1.7 

3.2 

2.8 

1.2 

2.0 

Sin(log) flattened T -only 

2.5 

0.7 

1.8 

2.5 

0.8 

1.9 

2.6 

0.9 

2.1 

3.0 

1.6 

3.2 

Sin(log) flattened T+E 

2.9 

1.4 

2.9 

3.0 

1.6 

3.4 

3.1 

1.6 

3.4 

3.6 

2.3 

4.5 


flattened resonance shapes offer stronger hints. The multi-peak 
equilateral signal rose from 1.9 cr (T-only) to 3.1 cr (T-\-E) af¬ 
ter adjusting for the “look-elsewhere” effect, while the flattened 
signal went from 2.4 cr (T-only) to 3.2 cr (T-\-E). These interest¬ 
ing results, reflecting those obtained for feature models, suggests 
the fit to any underlying NG signal might await alternative, but 
related, oscillatory models for a more compelling explanation. 

High frequency resonance model estimator: We have further sur¬ 
veyed the simple resonance model (Eq. 10) with a second ap¬ 
proach. Using a model-specific expansion in terms of linear os¬ 
cillation, proposed in Munchmeyer et al. (2015), it is possible to 
extend the frequency range of the analysis considerably. In this 
approach one exploits the fact that any bispectrum shape which 


is a function of ki+k 2 + h can be expanded in Fourier modes of 
ki -f ^2 + ^ 3 , resulting in an effectively one-dimensional expan¬ 
sion, as opposed to the general Modal expansion. In the present 
implementation we use 800 sine and cosine modes, which cover 
the full frequency space of the power spectrum search, i.e., 
0 < oj < 1100. The significances found with this method are 
presented in Fig. 25. As was the case for the feature model, we 
find that SMICA, SEVEM and NILC results are in good agreement 
in both temperature and polarization, while Commander results 
generally have larger variance and so are not included in the plot¬ 
ted averages. We find that the largest peak of the average is 3.6 cr 
in temperature and 3.1 cr in temperature and polarization com¬ 
bined. The results of the Modal expansion discussed in the pre- 
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vious paragraph are in good agreement with the high frequency 
resonance model estimator in the overlapping frequency range. 

Due to the high computational demands of this analysis, 
we have only exactly assessed the look-elsewhere effect for 
SMICA T data, for which we find an expected maximum peak 
of 3.5 (T + 0.4 (T in the case of Gaussian maps, to be compared 
to 3.7 (T in the SMICA T data, demonstrating that the results are 
fully consistent with Gaussianity. The expected maximum peak 
for Gaussian maps was calculated from the Fisher matrix with 
the method described in Meerburg et al. (2015). The average 
over component separation methods as well as the T+E data 
is even less significant. For the high frequency estimator we 
have assessed the significance of multiple peaks in the follow¬ 
ing way. Define the multi-peak amplitude Am as the sum of 
squares of the M highest significances ct/ in the frequency range, 

i.e.. Am = (Z^i ^ where only approximately indepen¬ 
dent peaks with Aoj > 10 are considered. One can then com¬ 
pare Am to its distribution in the Gaussian case and get an in¬ 
dividual significance ctm for each number of peaks M, where 
we assume M < 10. The multi-peak statistic is then obtained 
by maximizing over M, leading to an additional look-elsewhere 
effect that we also accounted for. In this way we find that in 
the SMICA T data the raw peak maximum 3.7 cr, is reduced 
to 0.5 (T for the single-peak statistic and to 0.6 a in the multi¬ 
peak case. This large reduction in significance is due to the 
large number of independent frequencies as well as to the max¬ 
imization over phases. One may argue that the frequency range 
oj < 1100 is too large, as EFT arguments for resonance non- 
Gaussianities (Behbahani et al. 2012) limit the frequency range 
to (9(10^). As an example, we have therefore also calculated the 
look-elsewhere-corrected significances when we limit the analy¬ 
sis to tu < 250. In this case we find a single-peak significance of 
0.6 (T and a multi-peak significance of 0.9 cr. Clearly the results 
are fully consistent with Gaussianity. 

8.4. Equilateral-type models and the effective field theory of 
inflation 

There is considerable interest in equilateral-type models because 
they are physically well-motivated, through e.g., varying sound 
speed scenarios. There are generic predictions available from 
the effective field theory of inflation, notably the two specific 
effective field theory (EFTl and EFT2) shapes that give rise 
to the equilateral and orthogonal approximations. These mod¬ 
els were previously constrained in Planck Collaboration XXIV 
(2014) and the reader is referred to section 2 of that paper for 
analytic expressions for the specific shapes constrained here. In 
Table 23, we list the main equilateral-type models in the liter¬ 
ature, giving constraints for T-only and T-\-E. All these models 
correlate well with the equilateral ansatz (Eq. 12) and likewise 
do not show a significant signal. However, despite this corre¬ 
lation, it is interesting to note the variation between models, 
largely due to the difference between these shapes in the flat¬ 
tened limit. The implications of these results are discussed in 
Sect. 11. 

8.5. Models with excited initial states (non-Bunch-Davies 
vacua) 

Non-Bunch-Davies (NBD) or excited initial states are mod¬ 
els which produce flattened (or squeezed) bispectrum shapes. 
The wide variety of NBD models that have been proposed are 
briefly classified and labelled in Sect. 2, following a more exten¬ 


sive overview in section 2 of Planck Collaboration XXIV (2014) 
where more analytic forms and the first constraints were pre¬ 
sented. The latest Planck constraints for these models are listed 
in Table 24, obtained using the Modal 2 estimator with polariza¬ 
tion, Despite the apparent “flattened” signal seen in the Planck 
bispectrum reconstructions (Fig. 4), this is generally not matched 
well by the specific modulation induced by the acoustic peaks 
for these scale-invariant NBD models. Tight constraints emerge 
for most models. The largest signal obtained is from the NBD 
sinusoidal shape which gives a 1.6 cr T-only raw significance, 
rising to 2.1 cr for T+E\ this is hardly an impressive correspon¬ 
dence given the number of models surveyed and the parameter 
freedom used in maximizing the signal. However, an important 
caveat for NBD models is that the predicted shapes can be very 
narrow in the flattened limit, in which case solutions have been 
smoothed to match the current Modal resolution (though this has 
improved considerably since the Planck 2013 NG analysis). An 
improved match to the warm inflation shape means that the final 
constraint shown in Table 24 is more robust, with further impli¬ 
cations discussed in Sect. 11 . 

8.6. Direction-dependent primordial non-Gaussianity 

We impose observational limits on direction-dependent primor¬ 
dial NG parametrized by Eq. 23. Rather than using ci and C 2 
we instead choose to work with the non-linearity parameters 
= -Cl/4 and = -C 2 /I 6 (chosen to match a primordial 
bispectrum that is equal to the equilateral shape in the equilateral 
limit) keeping the notation from the 2013 results. We estimated 
the values from temperature data and high-pass filtered po¬ 
larization data from the four foreground-cleaned CMB maps 
SMICA, NILC, SEVEM, and Commander, where we apply the com¬ 
mon mask. The details of the KSW estimator and its derivation 
is presented in Appendix A. For temperature data, we use the 
common mask as adopted in Planck Collaboration XII (2014), 
which has more conservative foreground masking than the newly 
available mask. We choose the more conservative foreground 
masking, considering the fact that anisotropic NG is more sen¬ 
sitive to residual foregrounds. We set the maximum multipole 
to 2000 and 1000 for temperature and polarization data, respec¬ 
tively. Validating our analysis pipeline with realistic simulations, 
we find that the asymmetry of the Planck beam, coupled with 
the Planck scanning pattern, inflates the statistical fluctuations 
of the significantly. Noting the large angular scale of arti¬ 
ficial anisotropy produced by the beam asymmetry, we set the 
minimum multipole to 101 , and find that the statistical fluctua¬ 
tion of estimation from simulated data is close to the theoretical 
expectations. 

These two shapes are also constrained using the Modal 2 es¬ 
timator, which is not affected by the beam asymmetry and is used 
in the same form as elsewhere in the paper with multipoles from 
2 to 2000 and 30 to 1500 being used for temperature and polar¬ 
ization, respectively. The present constraints are consistent with 
those found for T-only in Planck Collaboration XXIV (2014), 
but at higher resolution convergence has improved considerably, 
reflected in the lower variance. 

We find that the ISW-lensing bispectrum and the unresolved 
point-sources bispectrum bias the estimation of the in par¬ 
ticular in the analysis of temperature data. For our final val¬ 
ues, we subtract both these biases from our estimation. In Table 
25, we report the estimated value of from the foreground- 
cleaned CMB maps. For L = 1 the effect of the differing Aranges 
between the two estimators is not so significant and the results 
are quite consistent. For L = 2, which has significant signal in 
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Fig. 25. Standard resonance model results for both T-only and T+E across a wide frequency range using the high frequency resonance model 
estimator. The first and second panels show the signal in both T-only (first panel) and T+E (second panel) across the frequency range 0 < du < 500. 
The lower two panels give the results for T-only (third panel) and T+E (fourth panel) in the range 500 < io < 1100. 



the squeezed configuration, the effect of removing small scales 
from the KSW estimator is more pronounced, resulting in signif¬ 
icantly enlarged error bars. In light of this, the differences seen 
between the central values for the two methods is to be expected 
and does not indicate any inconsistencies between the two ap¬ 
proaches. The slight differences between the results from differ¬ 
ent foreground-cleaned temperature maps are within the likely 


range of statistical fluctuations, estimated from realistic simu¬ 
lations of CMB and noise propagated through the pipelines of 
foreground-cleaned map making. As seen in Table 25, we And 
that the estimated values of from Planck 2015 temperature 
plus polarization data are consistent with zero. 


47 




























Planck Collaboration: Planck 2015 Results. Constraints on primordial NG 

Table 23. Constraints on models with equilateral-type NG covering the shapes predicted by the effective field theory of inflation, 
together with constraints on specific non-canonical inflation models, such as DBI inflation. See Planck Collaboration XXIV (2014) 
(section 2) for further explanation of these specific models, with further implications discussed in Sect. 11. 



SMICA 


SEVEM 


NILC 


Commander 


Equilateral-type model 

A + cta 

S/N 

A + cta 

S/N 

A + cta 

S/N 

A + o-a 

S/N 

Constant T -only 

12 + 38 

0.3 

16 + 38 

0.4 

10 + 37 

0.3 

1 + 39 

0.0 

Constant T+E 

18 + 22 

0.8 

28 + 24 

1.2 

12 + 23 

0.5 

26 + 24 

1.1 

Equilateral T -only 

-15 + 68 

-0.2 

-9 + 68 

-0.1 

-19 + 67 

-0.3 

-17 + 69 

-0.3 

Equilateral T+E 

5 + 42 

0.1 

4 + 45 

0.1 

-2 + 42 

-0.1 

27 + 45 

0.6 

EFT shape 1 T-only 

-3 + 65 

0.0 

3 + 64 

0.0 

-7 + 62 

-0.1 

-9 + 66 

-0.1 

EFT shape 1 T+E 

12 + 39 

0.3 

15 + 42 

0.3 

3 + 39 

0.1 

32 + 42 

0.8 

EFT shape 2 T -only 

17 + 50 

0.3 

22 + 50 

0.5 

15 + 47 

0.3 

8 + 51 

0.2 

EFT shape 2 T+E 

23 + 29 

0.8 

31 + 31 

1.0 

15 + 29 

0.5 

36 + 31 

1.2 

DBI inflation T-only 

3 + 62 

0.0 

9 + 61 

0.1 

-1 + 59 

0.0 

-4 + 63 

-0.1 

DBI inflation T+E 

15 + 37 

0.4 

20 + 39 

0.5 

7 + 37 

0.2 

34 + 40 

0.9 

Ghost inflation T-only 

-50 + 80 

-0.6 

-45 + 80 

-0.6 

-54 + 79 

-0.7 

-45 + 82 

-0.6 

Ghost inflation T+E 

-27 + 50 

-0.5 

-37 + 54 

-0.7 

-31 + 51 

-0.6 

1 + 55 

0.0 

Inverse decay T -only 

17 + 43 

0.4 

21 +43 

0.5 

14 + 41 

0.3 

4 + 44 

0.1 

Inverse decay T+E 

23 + 25 

0.9 

32 + 27 

1.2 

15 + 26 

0.6 

32 + 27 

1.2 


Table 24. Constraints on models with excited initial states (non-Bunch-Davies models), as well as warm inflation. See Sect. 2 for 
further explanation and the labelling of these classes of NBD models. Note that the NBD, NBDl, and NBD2 models contain free 
parameters, so here we quote the maximum significance found over the available parameter range; the maximum for T and T+E can 
occur at different parameter values (on which the error bars are also dependent). 


SMICA SEVEM NILC Commander 


Flattened-type model 

A ± 

: CTa 

S/N 

A±(rA 

S/N 

A ± 

: CTa 

S/N 

A + 

CTa 

S/N 

Flat model T -only 

49 

+ 

65 

0.8 

57 + 

65 

0.9 

47 

+ 

65 

0.7 

19 

+ 

65 

0.3 

Flat model T+E 

44 

+ 

37 

1.2 

70 + 

37 

1.9 

33 

+ 

37 

0.9 

47 

+ 

37 

1.3 

Non-Bunch-Davies T-only 

42 

+ 

82 

0.5 

53 + 

82 

0.6 

26 

+ 

82 

0.3 

17 

+ 

82 

0.2 

Non-Bunch-Davies T+E 

61 

+ 

47 

1.3 

76 + 

47 

1.6 

43 

+ 

47 

0.9 

58 

+ 

47 

1.2 

NBD sine T-only 

567 

+ 

341 

1.7 

513 + 

341 

1.5 

588 

+ 

341 

1.7 

604 

+ 

341 

1.8 

NBD sine T+E 

-387 

+ 

206 

-1.9 

-485 + 

218 

-2.2 

-425 

+ 

206 

-2.1 

-417 

+ 

210 

-2.0 

NBDl cos flattened T-only 

-10 

+ 

22 

-0.5 

-4 + 

22 

-0.2 

-8 

+ 

22 

-0.4 

-9 

+ 

22 

-0.4 

NBDl cos flattened T+E 

-20 

+ 

19 

-1.1 

-10 + 

19 

-0.5 

-19 

+ 

19 

-1.0 

-14 

+ 

19 

-0.8 

NBD2 cos squeezed T-only 

10 

+ 

17 

0.6 

10 + 

17 

0.6 

8 

+ 

17 

0.5 

-2.5 

+ 

17 

-0.1 

NBD2 cos squeezed T+E 

-3 

+ 

5 

-0.5 

-0.8 + 

5.5 

-0.1 

-4 

+ 

5 

-0.8 

-3.8 

+ 

5.5 

-0.7 

NBDl sin flattened T-only 

-25 

+ 

22 

-1.1 

-27 + 

22 

-1.2 

-18 

+ 

22 

-0.8 

-33 

+ 

23 

-1.4 

NBDl sin flattened T+E 

48 

+ 

30 

1.6 

49 + 

33 

1.5 

35 

+ 

31 

1.1 

26 

+ 

34 

0.8 

NBD2 sin squeezed T-only 

-2.0 + 

1.4 

-1.4 

-1.4 + 

1.4 

-1.0 

-1.6 + 

1.4 

-1.1 

-1.3 

+ 

1.4 

-0.9 

NBD2 sin squeezed T+E 

-0.8 + 

0.4 

-1.9 

-0.5 + 

0.4 

-1.2 

-0.6 + 

0.4 

-1.4 

-0.5 

+ 

0.4 

-1.2 

NBD3 non-canonical T-only (xlOOO) 

-5.9 + 

6.7 

-0.9 

-6.0 + 

6.8 

-0.9 

-5.4 + 

6.8 

-0.8 

-5.5 

+ 

6.7 

-0.8 

NBD3 non-canonical T+E (xlOOO) 

-8.7 + 

5.0 

-1.7 

—6.2 + 

5.2 

-1.2 

-7.5 + 

5.2 

-1.5 

-9.4 

+ 

5.2 

-1.8 

Warms inflation T -only 

-23 

+ 

36 

-0.6 

-26 + 

36 

-0.7 

-32 

+ 

36 

-0.9 

-24 

+ 

36 

-0.7 

WarmS inflation T+E 

-14 

+ 

23 

-0.6 

-28 + 

23 

-1.2 

-21 

+ 

23 

-0.9 

-17 

+ 

23 

-0.7 


8.7. Parity-violating tensor non-Gaussianity 

We present observational limits on the parity-violating tensor 
nonlinearity parameter from the temperature and £’-mode 
polarization data. Unlike the usual scalar-mode templates, the 
CMB bispectra sourced from the tensor NG (Eq. 27) are written 
in non-factorizable forms (Shiraishi et al. 2013b); hence, we use 
the separable Modal pipeline in our bispectrum estimations. 

The parity-violating NG under examination induces non¬ 
vanishing signals not only in parity-even configurations (^i -r 
^2 + ^3 = even) but also in the parity-odd ones (^i -r ^2 + 
i'i = odd) in the temperature and £’-mode polarization bispec¬ 
tra (Shiraishi et al. 2013b). The optimal estimator, including all 
(even -r odd) bispectrum signals, is expressed by the linear com¬ 
bination of the parity-even and parity-odd estimators, reading 


(Liguori et al. 2015) 

yw^even ^even , Ajodd fodd 

/all = /nL 

/yeven _|_ /yodd ^ 

where A/^^ven/odd normalization factor (related to the 

Fisher matrix as A^even/odd _ ^^even/odd^ defined for T’l -h ^’2 + 
^3 = even/odd. The parity-even estimator can be com¬ 
puted using the original Modal methodology (Fergusson et al. 
2010a, 2012; Fergusson 2014; Figuori et al. 2015), while in 
computations of the parity-odd estimator we follow the 
extended spin-weighted pipeline (Shiraishi et al. 2014, 2015; 
Figuori et al. 2015). 

Our estimations (with both temperature and polariza¬ 
tion data) are based on the resolution of ^max = 500 and HEALPix 
A^side = 512, leading to feasible computational costs. This choice 
is not expected to change the results significantly, in compari- 
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Table 25. Direction-dependent NG results for both the L = \ and L = 2 models. We present results from both the KSW and Modal 
2 pipelines. The discrepancy between the central values for the L = 2 models is due to the differing i ranges taken for the two 
estimators, the key difference being the KSW ^min = 101. As this model peaks in the squeezed configuration, a significant portion 
of the signal is lost, which is reflected in the increased error bars. 



Commander 


NILC 


SEVEN 


SMICA 



A±(Ta 

S/N 

A + cta 

S/N 

A + cta 

S/N 

A±(Ta 

S/N 

L= 1 

Modal 2 T -only 

-41 + 43 

-0.9 

-58 + 42 

-1.4 

-51 +43 

-1.2 

-49 + 43 

-1.1 

KSW T-only 

-8 +46 

-0.2 

-62 + 46 

-1.3 

-34 + 45 

-0.8 

-26 + 45 

-0.6 

Modal 2 T+E 

-28 + 29 

-1.0 

-30 + 27 

-1.1 

-49 + 28 

-1.7 

-31 +26 

-1.2 

KSW T+E 

-57 + 33 

-1.7 

-62 + 32 

-1.9 

-79 + 32 

-2.5 

-54 + 32 

-1.7 

L = 2 

Modal 2 T -only 

0.7+ 2.8 

0.2 

0.8+ 2.8 

0.4 

1.1 + 2.7 

0.3 

0.5+ 2.7 

0.2 

KSW T-only 

1.5+ 5.1 

0.3 

-3.9+ 5.1 

-0.8 

-0.4+ 5.1 

-0.1 

0.1+ 5.0 

0.0 

Modal 2 T+E 

1.1 + 2.4 

0.5 

0.5+ 2.4 

0.2 

1.3+ 2.4 

0.6 

-0.2+ 2.3 

-0.1 

KSW T+E 

-3.0+ 4.1 

-0.7 

-3.6+ 4.0 

-0.9 

-3.8+ 4.0 

-1.0 

-1.3+ 3.9 

-0.3 


son to the analysis at higher resolution, e.g., ^max = 2000 and 
Healpix Aside = 2048, since the cosmic variance and instru¬ 
mental noise are already far higher than the signals for ^ > 300 
(Shiraishi et al. 2013b). Only in the polarization data analysis is 
an effective ^min also adopted, which is motivated by the high- 
pass filtering process for ^ < 40 in component separation. 

Within the above I ranges, the theoretical bispectrum tem¬ 
plates are decomposed with the eigenbasis composed of (9(1 - 
10 ) polynomials and some special functions reconstructing the 
tensor-mode features, e.g., temperature enhancement due to the 
ISW effect {i < 100 ), and two £’-mode peaks created by reion¬ 
ization {i < 10) and recombination {i ^ 100). The resulting 
factorized templates are more than 95 % correlated with the orig¬ 
inal ones. The validity of our numerical computations has been 
confirmed through the map-by-map comparisons of at 

very low resolution, showing the consistency between the values 
from the Modal methodology and those obtained by the brute- 
force 0{l^) summations like Eq. (36). We have also checked that 
our parity-even estimator successfully leads to the constraints on 
^max = 500, Compatible with the results 
from the binned estimator. 

Our limits estimated from the foreground-cleaned tempera¬ 
ture and high-pass filtered polarization data (SMICA, SEVEN, and 
NILC) are summarized in Table 26. The data and MC simula¬ 
tions used here, including all experimental features, i.e., beam, 
anisotropic noise levels and partial sky mask, have been in- 
painted using the identical recursive process adopted in the stan¬ 
dard /nl estimations (see Sect. 3.5). The sky fractions of the tem¬ 
perature and polarization masks adopted here are, respectively, 
/sky = 0-26 and /sky = 0.74. Although the error bars and the 
linear terms have been computed using 160 MC simulations, the 
resulting error bars are close to the expected values, (/skyE)“^/^. 

We have confirmed the stability of the T-only constraints, 
and significant scatter of the £’-only constraints both in the 
parity-even case and in the parity-odd one, when changing /ky. 
Such £’-mode instability has given insignificant effects on our 
T+E constraints in the parity-even case, as they are determined 
almost exclusively by TTT, like the scalar NG analyses. In con¬ 
trast, our parity-odd T+E results vary a lot, due to the £'-mode 
scatter (quantitatively speaking, only a few percent change of 
/sky has shifted by more than 1 cr), because TTE and TEE 
contribute significantly to the signal-to-noise ratio in the parity- 
odd case (Shiraishi et al. 2013b). Table 27 presents the corre¬ 
lations of the bispectra reconstructed from the component sep¬ 


arated maps, also indicating the robustness of the T-only con¬ 
straints and the instability of the £’-only results. We report only 
stable results in Table 26 and conclude that there is no evidence 
at > 2 cr of in the parity-even, parity-odd or their whole 
domain. 

The parity-odd components of the TTT and EEE bispectra 
extracted model-independently from the SMICA data are visually 
represented in Fig. 26. It is apparent from this figure that the 
Planck TTT bispectrum has similar features to the WMAP one 
(Shiraishi et al. 2015), e.g., distinctive signals distributed around 
£i ^ ~ ^ 3 . As indicated by the roughly 70% correlation 

between the SMICA and WMAP bispectra (see Table 27), the 
Planck T-only limits in Table 26 are close to the WMAP ones 
(68 % CL): = 4+16 for parity-even; and = 

80+110 for parity-odd (Shiraishi et al. 2015). 

Table 26. Results for the tensor nonlinearity parameter 
estimated from the SMICA, SEVEN, and NILC temper¬ 
ature and high-pass filtered polarization maps. We separately 
show the central values and the errors (68 % CL) extracted from 
ii + £2 + h = even (Even), + ^2 + ^3 = odd (Odd) and their 
whole domain (All). The parity-odd constraints have also been 
obtained from the £’-mode data, but they are still preliminary and 
not currently shown. 


Even Odd All 


SMICA 

T . 2+ 15 120+ 110 4+ 15 

T+E . 0+13 

SEVEN 

T . 2+ 15 120+ 110 5 + 15 

T+E . 4+13 

NILC 

T . 3 + 15 110+ 100 5 + 15 

T+E . 1 + 13 


9. Limits on the primordiai trispectrum 

So far, we have considered a variety of physically motivated 
possibilities for the inflationary 3-point function, or bispectrum. 
A similar phenomenology exists for the 4-point function, or 
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trispectrum. Our constraints on the trispectrum will use CMB 
temperature only; we do not use polarization in this section. We 
start by briefly reviewing the inflationary physics and classifying 
the signals we will search for. 

First, some notation: a “primed” f-trispectrum 
{^ki^k 2 ik 3 ^kAy denotes the connected trispectrum without 
its momentum-conserving delta function, i.e., 

(^ki(k 2 (kA) = 2 ki) + disc., (62) 

where “-r disc.” denotes the disconnected contributions to the 
4-point function. 

One possible signal is the “local” trispectrum which 
arises if the non-Gaussian adiabatic curvature f is of cubic-type 
form (see, e.g., Okamoto & Hu (2002)), i.e., 

^ix) = ^g(-^) + ^ (63) 



where is a free parameter and is a Gaussian held. In this 
model, the bispectrum is zero (since there is a f ^ symme¬ 
try) and the trispectrum is given by 

{^k.^k^^kAY = § 0Nr‘ [Pdh)P((k2)P((k3) + 3 perm.]. (64) 

Analogously to the case of the local bispectrum, 
the observational signal-to-noise for g^^^^ is largest in the 
“squeezed” limit, ki «: min(^ 2 ,^ 3 ,^ 4 ), and there is a con¬ 
sistency relation which shows that in single-fleld inflation, 
the four-point function is always small in the squeezed 
limit (e.g., Senatore & Zaldarriaga 2012a). Thus g^^^^ can only 
be detectably large in multi-fleld models. Conversely, there are 
multi-fleld models where g^^^^ is detectable. The main obstacle 
here is technical naturalness, i.e., ensuring that radiative correc¬ 
tions do not generate an observationally larger bispectrum. This 
can be the case if a large bispectrum is forbidden by a Z 2 sym¬ 
metry, or by supersymmetry (Senatore & Zaldarriaga 2012b). 

A further category of four-point signals can be obtained 
by adding quartic interactions to the inflationary action. 
Following Smith et al. (2015), we will concentrate on the sim¬ 
plest possibility, by considering quartic operators consistent with 
the symmetries of inflation and having the lowest possible num¬ 
ber of derivatives (Bartolo et al. 2010b; Senatore & Zaldarriaga 
2011; Senatore & Zaldarriaga 2012b). There are three such oper¬ 
ators, of the schematic form (t"^, &^(dicr)^, and (dicr)^(djcr)^. By 
a short calculation using the in-in formalism (Maldacena 2003), 



Fig. 26. Parity-odd signals (^1 -r ^2 + ^3 = odd) of the TTT 
(top) and EEE (bottom) bispectra {it < 500) recovered from 
the SMICA maps by means of the Modal decomposition with 101 
simple polynomial-based eigenmodes, not including any special 
functions fitting the CMB tensor-mode features. In the panel for 
EEE, only the signals larger than ^ = 40 are shown. The TTT and 
EEE bispectra shown here are rescaled with a constant Sachs- 
Wolfe weighting and signal-to-noise weighting, respectively. 


Table 27. Correlation coefficients between pairs of bispectrum 
modes, extracted from two of the Planck component separated 
maps and the WMAP foreground-cleaned map at ^max = 500 
resolution. We separately present the results estimated from -r 
^2 + ^3 = even (Even) and ^1 -r ^2 + ^3 = odd (Odd) combinations. 
The loss of the correlations is confirmed in the EEE case, like 
Table 14. 



TTT 

EEE 

Methods 

Even 

Odd 

Even 

Odd 

SMICA-SEVEM . . 

1.00 

0.99 

0.80 

0.80 

SMICA-NILC . .. 

1.00 

1.00 

0.90 

0.87 

SEVEM-NILC . . . 

0.99 

1.00 

0.70 

0.60 

SMICA-WMAP . 

0.75 

0.67 




the associated four-point functions are (Smith et al. 2015; see 
also Huang & Shiu 2006): 




25 
221184 
25 


g^^A] I dTE^ n 




^NL 


1 


^1^2 ^3 ^4^^ ’ 


(65) 
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// / / r - 13824 or^ido-f . 3 

\hk\hk2bk^hkA/ — ^NL 


325 


drsTl 


(1 -ksTsXl -k4TE) 


27648 


k\k2k^k^ 
+5 perm 

,o-\do-f .3 


(*3 • k4)e^ 


kiTE 


325 

+ 3(^3 + ^4)^ + 12^3^4 


X 


^1^2 ^3 ^4^^ 

+5 perm.j; 


(A:3 • k^) 


( 66 ) 


\bkihk2bk^hk^) — ntcnc ^NL 


2575 


f 


cIte 


n 

1=1 


(1 - fc, T£)e‘ 


kiTE 




X [(A:i • A:2)(A:3 • A:4) + 2 perm.] 

165888 (5fr)%3 
2575 

2/:4 -2K^Y.kJ + KY,k^^ + \2k1k2k3k4 


k\k^ky<?^K^ 

X [(A:i • A:2)(A:3 • A:4) + 2 perm.] 


(67) 


Here K = Yjh and we have introduced parameters ’ 

and to parametrize the amplitude of each trispectrum. 

The normalization of the ^NL-parameters is chosen so that 
{^kA 2 ^k,^k,) = (216/25) for “tetrahedral” configura¬ 

tions, with \ki\ = k and {ki • kj) = -k^/3. This is the analogue 
of the commonly-used normalization for the bispectrum, where 
/nl parameters are defined so that {^ki^k 2 ^k 3 ) = (18/5) 
for equilateral configurations with \ki\ = k. 

For simplicity in Eqs. (65)-(67) we have assumed a scale- 
invariant initial power spectrum P^ik) = A^jk^. In order to anal¬ 
yse Planck data, we must slightly generalize this to a power- 
law spectrum P^(k) oc Our scheme for doing this follows 
Appendix C of Smith et al. (2015). 

A Fisher matrix analysis shows that there is one large corre¬ 
lation among the three trispectra in Eqs. (65)-(67), so that to an 
excellent approximation we can treat only two of the trispectra as 
independent. To quantify this, in Smith et al. (2015) it is shown 
that the &^(dcr)^ shape is 98.6 % correlated to a linear combina¬ 
tion of the shapes and (dcr)'^. Therefore, we will only search 

for the parameters and g^^[^ . 

We note that the analysis that leads to the trispectrum shapes 

is similar to the analysis that leads to the 

“standard” bispectrum shapes and However, there 

are some minor differences as follows. In the bispectrum case, 
one considers the cubic operators and 7i{d7if, but it is con¬ 
ventional to define observables which are related to 

the operator coefficients by a linear transformation. This is done 
because the two cubic operators are about 90 % correlated, so it 
is convenient to orthogonalize. In the trispectrum case the corre¬ 
lation is smaller (around 60 % for Planck), and we have chosen 
to omit the orthogonalization step. Another reason to omit the or- 
thogonalization step is that the trispectrum shape (dcr)'^ is a sig¬ 
nature of multi-field inflation. In single field inflation, the (dcr)^ 


trispectrum is not technically natural; radiative corrections gen¬ 
erate cubic operators of the form or nidn)^, which generate a 
bispectrum with larger signal-to-noise. 

There are more trispectrum shapes one might consider. For 
example, classifying Galilean invariant quartic operators leads 
to higher-derivative trispectra, which are not highly correlated to 
the trispectra considered above (Bartolo et al. 2013; Arroja et al. 
2013). We have only considered “contact” diagrams arising from 
one power of a quartic operator, and it would be interesting 
to study “exchange” diagrams arising from two cubic opera¬ 
tors and exchange of a light particle (e.g., Chen et al. 2009; 
Arroja et al. 2009; Chen & Wang 2010; Bartolo et al. 2010a; 
Baumann & Green 2012). We leave these as extensions for fu¬ 
ture work. 

Summarizing, we will search for the following trispectrum 
signals: 

(68) 

defined by Eqs. (64), (65), and (67) above. 

9.T Dafa analysis 

Turning now to data analysis, we use the machinery 
from Smith et al. (2015). The first step is to represent each 
trispectrum as a small sum of factorizable terms as follows. The 
angular CMB trispectrum can be written either as an integral 
over comoving distance r (in the case of or as a dou¬ 

ble integral over (r, r) where r is conformal time during infla¬ 
tion (in the case of the or (dcr)"^ trispectra). We approximate 
the integral by a finite sum, which represents the CMB trispec¬ 
trum as a sum of terms that are factorizable in a sense defined 
in Smith et al. (2015). A large number of sampling points are 
needed to obtain a good approximation to the integral, lead¬ 
ing to a large number of terms in the factorizable representa¬ 
tion. However, there exists an optimization algorithm, which 
takes as input a trispectrum that has been represented as a sum 
of many factorizable terms, and outputs a representation with 
fewer terms. The reduction can be quite dramatic, as shown in 
Table 28. The optimization algorithm guarantees that the output 
trispectrum accurately approximates the input trispectrum, in the 
sense that the two are nearly observationally indistinguishable. 


Table 28. Number of factorizable terms Ain needed to represent 
each trispectrum by direct sampling of the integral, and number 
of terms Aout obtained after running the optimization algorithm 
from section VII of Smith et al. (2015). 


Trispectrum 


Aout 

„local 

^NL . 

436 

17 

. 

. 6955 

73 

{dcrf . 

. 20865 

192 


Armed with “small” factorizable representations for each 
trispectrum, the next step is to run an analysis pipeline that 
estimates the amplitude of each trispectrum from Planck data. 
We use the “pure MC” pipeline from Smith et al. (2015), which 
compares the trispectrum of the data to the mean trispectrum 
of an ensemble of simulations. This pipeline requires a filtering 
operation d dim which processes the pixel-space CMB data d 
and generates a harmonic-space map dim- Otir filtering operation 
is defined by the following steps: 
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1. Starting from the data d, we compute (with uniform 
pixel weighting) a best-fit monopole and dipole outside 
the Galactic mask. We use the temperature “common 
mask”, the union of the confidence masks for the SMICA, 
SEVEN, NILC, and Commander component separation meth¬ 
ods (Planck Collaboration IX 2016). 

2. The mask defines a few “islands”, i.e., isolated groups of 
pixels that are unmasked, but contained in a larger masked 
region. We slightly enlarge the mask so that it removes the 
islands. 

3. We classify the components of the masked part of the sky 
into “small” masked regions with < 1000 pixels (at HEALPix 
resolution A^side = 2048), and “large” regions with > 1000 
pixels. Small regions usually correspond to point sources, 
and large regions typically correspond to areas of diffuse 
galactic emission. In small regions, we inpaint the CMB by 
assigning the unique map that agrees with the data on bound¬ 
ary pixels, and whose value in each interior pixel is the aver¬ 
age of the neighboring pixels. 

4. In large regions, we do not inpaint the CMB, but rather 
apodize the boundary of the large region using cosine 
apodization with 12’ radius. 

5. We apply a spherical harmonic transform to the inpainted, 
apodized CMB map to obtain a harmonic-space map 
with ^max = 1600. We then take the final filtered map 

to be 


— 


b(Cc + b~^Nc 


(69) 


where is the beam, Q is the fiducial CMB power spec¬ 
trum, and Ni is the sky-averaged noise power spectrum 
(without beam deconvolution). To motivate this choice of i- 
weighting, we note that for an ideal all-sky experiment with 
isotropic noise, we have a^rn = where are 

signal and noise realizations. In this case, Eq. (69) weights 
the signal as ^^m/(Q + which is optimal. 


In our pipeline, we apply this filter to the component- 
separated SMICA maps (Planck Collaboration IX 2016), obtain¬ 
ing a harmonic-space map We apply the same filter to 1000 
Monte Carlo simulations to obtain an ensemble of harmonic- 
space maps. Our pipeline has the property that it always esti¬ 
mates the trispectrum of the data in excess of the trispectrum 
in the simulations. Since the simulations include lensing, this 
means that lensing bias will automatically be subtracted from 
our ^nl estimates. 

Now that the filter, data realization, and Monte Carlo simula¬ 
tions have been fully specified, the details of the pipeline are de¬ 
scribed in section IX.B of Smith et al. (2015). For each trispec¬ 
trum, the pipeline outputs an estimate of ^nl and an estimate of 
the statistical error. Our basic results are: 


Our central value in Eq. (70) agrees well with this result, but the 
statistical error is smaller by a factor of 2.3. This improvement 
is partly due to the lower noise levels in Planck 2015 data, and 
partly due to the use of a better estimator. 

Each line in Eq. (70) is a “single-^NL” constraint; i.e., the 
constraint on one ^nl parameter with the other ^nl- parameters 
held fixed. For joint constraints, one needs to know the full co- 
variance matrix. The correlation between and the other two 

parameters is negligble, and the correlation is: 

Corr(g^,g^^[^") = 0.6l. (71) 

multi-field models of inflation will generally give a linear com¬ 
bination of (T^, &^(dicr)^, and (di(T)^(dj(T)^ trispectra. In this case 
we proceed as follows. First, if the &^(dicr)^ coefficient is non¬ 
zero, we can use the near-degeneracy with a linear combination 
of the other two operators to absorb it into the effective values 

of ^nl ^ ^ Fisher matrix analysis shows that the coeffi¬ 

cients of this linear combination are 

(<)eff = 

= 0.091 (72) 

It is convenient to define the two-component parameter vector: 

„=(<4 

We also compute a two-by-two Fisher matrix Fij, whose diago¬ 
nal is given by Fa = l/cr^, where ct/ is the single-^NL statistical 
error in Eq. (70), and whose off-diagonal is F 12 = 
where r is the correlation in Eq. (71). This procedure gives: 


FiJ = 


3.3 9.2 
9.2 68.7 


X 10 


■13 


(74) 


For a given parameter vector we can define a trispectrum-;p^ 
by 

x\g) = [Fugi - (Fg)i] F-/ [Fjjgj - (Fg)j] (75) 

where gt = (-0.21 x 10^,-0.10 x 10^) is the vector of best-fit 
single-^NL values from Eq. (70). This definition of follows 
from the observation that (Fagi) is an estimator with expectation 
value (Fg)i and covariance matrix Cow^Fugi, Fjjgj) = Fij. 

The inflationary implications of these trispectrum constraints 
are discussed in Sect. 11.5 below. 


10. Minkowski functionais resuits 


= (-9.0 ± 7.7) X 10“^; 

= (-0.2±1.7)xl0®; (70) 

9^^' = (-0.1±3.8)xlO^ 

No deviation from Gaussian statistics is seen. These results 
significantly improve the previous best constraints on the 
trispectrum from WMAP (Vielva & Sanz 2010; Smidt et al. 
2010; Fergusson et al. 2010b; Hikage & Matsubara 2012; 
Sekiguchi & Sugiyama 2013; Regan et al. 2013; Smith etal. 
2015) and large-scale structure (Desjacques & Seljak 2010; 
Giannantonio et al. 2014; Leistedt et al. 2014). 

A constraint on g^^^^ from Planck 2013 data was recently 
reported by Feng et al. (2015), who find = (-13 +18)x 10"^. 


In this section, we present constraints on local NG at 
first and second order and obtained with 

Minkowski functionals (MFs) on temperature and polarization 
E maps. MFs (Mecke et al. 1994; Schmalzing & Buchert 1997; 
Schmalzing & Gorski 1998; Winitzki & Kosowsky 1998) are a 
measure of fields’ local morphology used to constrain their 
stationarity, isotropy and Gaussianity. Mostly probing general 
NG in a frequentist fashion in two-dimensions on CMB maps 
(Eriksen et al. 2004; Komatsu etal. 2005; Modest etal. 2013; 
Natoli et al. 2010; Curto et al. 2008) or three-dimensions on 
ESS data (Park etal. 2005; Wiegand et al. 2014), they have 
also been used to measure specific NG targets with Bayesian 
methods, such as (Hikage et al. 2006, 2008; Ducout et al. 
2013; Planck Collaboration XXIV 2014), other bispectrum and 
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trispectrum shapes (Hikage & Matsubara 2012) and topologi¬ 
cal defects (Planck Collaboration XXV 2014). New develop¬ 
ments have been made recently, using needlets (Fantaye et al. 
2015), neural networks (Novaes et al. 2014) or allowing scale- 
dependent measurements (Munshi et al. 2013). 

MF-based limits are well known to be suboptimal for 
and but they provide an independent cross-check of bis¬ 
pectrum and trispectrum-based estimators. They are comple¬ 
mentary to optimal estimators: they are weighted integrals of 
the poly spectra and are sensitive to any source of NG, including 
foregrounds and secondaries. While MFs are not always able to 
distinguish between these different sources and systematics, they 
allow upper bounds to be put on them. 

The most recent constraints (1 cr) from MFs on and 
have been obtained with WMAP (Hikage & Matsubara 
2012) and Planck (Planck Collaboration XXIV 2014): 


From these measurements, we then use a Bayesian method 
to jointly estimate and 


p//) I flocal JocaKp/flocal Jocal\ 
r I fiscal JocaKp/flocal JocaKp flocal^Jocal ’ ^ '' 

j i7nl ’ ^NL ’ ^NL 

We take a uniform prior for in the range -400 to 400, and 
for g^^^^ in the range -4 x 10^ to -1-4 x 10^, while the evidence is 
just considered as a normalization. 

Assuming MFs are multi-variate Gaussian-distributed we 
obtain the posterior distribution for with test 




Jocal JocaK 

A yy^J^L ’^NL ' 


(80) 


=4.2 + 20.5-, = (1.9 + 6.4)x 10^ (76) 

10.1. Method and definition of MFs 

For a smoothed two-dimensional field S of zero mean and of 
variance cTq, defined on the sphere, we consider an excursion 
set of height v = S/ctq, i.e., the set of points where the field 
exceeds the threshold v. We use four functionals denoted by 
Vk(y)(k = 0,1,2,3). The first three correspond to MFs: Vo is 
the fractional Area of the regions above the threshold, Vi is the 
Perimeter of these regions and V 2 is the Genus, defined as the to¬ 
tal number of connected components of the excursion above the 
threshold minus the total number of connected components un¬ 
der the threshold. The fourth, V 3 , is the Number of clusters (also 
referred to as Ndusters)- This is the number of connected regions 
above the threshold for positive thresholds and below the thresh¬ 
old for negative thresholds. Precise definitions and formulae for 
the quantities Vk as well as their expectation values for Gaussian 
fields are summarized in Appendix B. 

We calculate the four normalizedfunctionals Vk(y) on Uth = 
26 thresholds v, between Vmin = -3.5 and Vmax = +3.5. 

For this analysis, we used the same temperature and polar¬ 
ization E data, simulations and masks described in Sect. 3.4 
for consistency with the bispectrum estimators. In addition, 
the maps are filtered to optimize constraints on local NG 
(Ducout et al. 2013), the filters used being similar to Wiener fil¬ 
ters for T and E (Wm), and for the first (Wdi) and second (IVd 2 ) 
derivatives of these fields (Fig. 27): 

Wdi k ^j({(+\)Wu■, (77) 

Wd 2 oc /(/+1)1%. (78) 

For the temperature map, known point sources in the mask are 
inpainted. 

We define the vector g as any combination y = with 

k = {0,3}, A = [T, E), W = {Wm, Wm, Wm}, y being the vector 
measured on the data. 


Raw Minkowski functionals Vk depend on the Gaussian part of 
fields through a normalization factor Ak, a function only of the shape 
of the power spectrum. We therefore normalize functionals Vk = Vk/^k 
to focus on NG, see Appendix B. 


with 

“ ^siml(/NL^\ ^NL^^)] ^sim2 “ ^siml(/NL^\ ^NL^^)] • (^1) 

For this test, we use two types of simulations to first con¬ 
struct a model including primordial NG ^simi(/NL^\ and 
secondly a covariance matrix 

Gsim2 = ^(?/sim2 “ ^sim2) (?/sim2 “ ^sim2) ^ ? (^2) 

with ^sim = (?/sim)sim averaged over the simulations. We now 
describe the details of these simulations. 

- Simulations 1: Non-Gaussian model 

For the first type of simulation, we included all possible 
sources of NG, assuming that the total and individual lev¬ 
els of NG are small enough that MFs are linear with respect 
to those NG levels (Ducout et al. 2013). The three kinds of 
NG we included are foreground residuals (Galactic residu¬ 
als with scalable amplitude a, as well as radio sources, GIB 
anisotropies, secondaries (SZ, lensing and ISW-lensing, but 
not SZ-lensing) and primordial NG (/nl^\ 

siml' = 

+ map^^(radio sources, GIB, SZ) 

+ a X map^^(Galactic residuals). (83) 

We tried to reproduce all instrumental effects, with realis¬ 
tic effective beams (isotropic window functions), noise from 
FFP8 simulations (Planck Collaboration XII 2016), filtered 
with component separation weights. We checked the ac¬ 
curacy of these simulations by comparing them to FFP8 
simulations, using no foreground and no primordial NG. 
The astrophysical models are provided by the Planck Sky 
Model (PSM, Delabrouille et al. 2013), while the primor¬ 
dial NG simulations are computed as in Eisner & Wandelt 
(2009). The lensing uses LensPix^^. The power spec¬ 
trum used for these NG simulations and the lensing is the 
best-fit power spectrum from Planck 2013+ACT/SPT+BAO 
(Planck Collaboration XXII 2014). We created ni = 200 
simulations i, using ni maps for the primordial NG, while 
we had only one astrophysical foreground simulation. 

http://cosmologist.info/lenspix 
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Fig. 27. Filters used to optimize constraints on local NG, in har¬ 
monic space. The temperature filter Wu is a smoothed version 
of the true Wiener filter obtained with realistic models, while the 
£’-mode Wu filter is adapted from the temperature one, with a 
cutoff value at ^ ^ 800. The formulae for the derivative filters 
are given in Eq. (78). 


- Simulations 2: FFP8 (Planck Collaboration XII 2016) MC 
simulations 

Since NG is weak, the covariance matrix C is computed with 
n 2 = 10"^ simulations, including no primordial NG and no 
foregrounds. These simulations reproduce realistic instru¬ 
mental effects (anisotropy of beams in particular), realistic 
noise and component separation filtering. The only NG still 
present in these simulations are lensing and the ISW-lensing 
correlation. 


Validation of the estimator 

Part of the validation for the MFs estimator is described in 
Sect. 5.3 for to compare the results to bispectrum estima¬ 
tors on realistic simulations (FFP8 MC, second item above). In 
addition we present in Table 29 the results obtained on the same 
realistic simulations for and on simulations containing pri¬ 
mordial NG (first item above), with = 10 and = 10^, 
with 200 simulations used in each case. These tests have been 
performed using the SMICA method with lensing bias removed. 

10.2. Results 

Results for and g^^^^ estimation with MFs on the four com¬ 
ponent separated maps in temperature and polarization are pre¬ 
sented in Table 30. The results for polarization £’-only maps are 
not quoted, since these results were not sufficiently stable (cf. 
Sect. 7.6). No deviation from Gaussianity is detected. T+E anal¬ 
ysis generally finds higher values for but remains consis¬ 
tent with Gaussianity. 

The posteriors for and g^^^^ from SMICA are shown in 
Fig. 28. One interesting point is that the estimates of 
g^^^^ are almost uncorrelated (r < 0.1); this can be inferred when 
we consider the parity of MF deviations from Gaussianity, which 
is different for the two parameters (Matsubara 2010). 

Foreground and secondary biases are removed from these es¬ 
timates, since the NG model directly includes them. However, 
an estimation of their contribution in the map is reported in 
Table 31. 


Table 29. Results for local NG parameters at first and second 
order, and obtained with Minkowski functionals on 
SMICA simulations in temperature and polarization. These re¬ 
sults are corrected for the lensing and ISW-lensing biases unless 
stated otherwise. Parameters are estimated jointly, and we report 
marginalized results, quoting 1 cr errors. The results are the av¬ 
erage obtained from 200 simulations. 




/•local 

-/NL 

(xio-*) 

/•local _ n „local 

JnL ~ ^ ’ ^NL 

T . 

= 0 

0+13 

-1 + 19 

E . 


1 +42 

0 + 23 

T + E . 


1 + 12 

0+13 



/•local — 1 A ^local 
JnL ~ ’ ^NL 

T . 

= 10^ 

10 + 13 

9 + 22 

E . 


12 + 42 

10 + 23 

T + E . 


11 + 12 

10+13 




Table 30. Results for local NG parameters at first and second 
order, and obtained with Minkowski functionals on 
all four component separated maps in temperature and polariza¬ 
tion. These results are corrected for the lensing and ISW-lensing 
biases unless stated otherwise. Parameters are estimated jointly, 
and we report marginalized results, quoting 1 cr errors. 



/•local 

-/NL 

(xlO") 

SMICA 

T . 

2+13 

-17 + 19 

T + E . 

3 + 12 

-8 + 13 

SEVEM 

T . 

3 + 13 

-23 + 20 

T + E . 

7 + 12 

-9 + 13 

NILC 

T . 

10+13 

-23 + 20 

T + E . 

12 + 12 

-15 + 13 

Commander 

T . 

8 + 13 

-30 + 19 

T + E . 

10+13 

-18 + 13 


Foreground and secondary biases 

Foreground residuals are generally negligible, in particular in the 
T analysis. This is different from the Planck 2013 results, where 
the residuals were more important; this can be explained by the 
beam correction applied to these previous estimates which exag¬ 
gerated signals from small scales. 

Fensing has a significant signature in MF estimation of 
but is even stronger in g^^^^ (the four-point correlation signature) 
and could be detected (and not treated just as a bias) with this 
estimator. The Wiener filters enhance the scales where lensing is 
dominant. 
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Fig. 28. Joint constraint on and obtained with MFs. 
The contour lines represent 1, 2 and 3 cr limits of a 2D-Gaussian 
distribution. Constraints were obtained with SMICA temperature 
and polarization E maps. 


Table 31. Biases for local NG parameters at first and sec¬ 
ond order obtained with Minkowski Functionals on 

SMICA in temperature and polarization. Parameters are estimated 
jointly, and we report marginalized results. For the correspond¬ 
ing error (on one map), see Table 30. 


Ag}^^i(xl04) 


T T+E T T+E 


SMICA 


SZ. 0.0 -0.3 2.3 1.1 

CIB . 0.7 0.5 -6.8 3.4 

Galaxy . -0.1 -0.2 -1.2 3.1 

PS . 0.1 0.2 2.2 1.2 

Lensing . 16.5 10.0 63.1 40.4 


11. Implications for eariy Universe physics 

The NG constraints obtained in this paper show consis¬ 
tency of Planck data with Gaussian primordial fiuctua- 
tions, thus confirming the results obtained in the 2013 re¬ 
lease (Planck Collaboration XXIV 2014) and improving them 
through the inclusion of CMB polarization data. The standard 
single-field slow-roll models of infiation have therefore been 
confirmed as a viable scenario for infiation, passing one of their 
most stringent tests, based on lack of measurable deviations 
from Gaussianity. The constraints obtained on local, equilat¬ 
eral, and orthogonal NG, after accounting for various contam¬ 
inants, strongly limit different mechanisms proposed as alterna¬ 
tives to the standard inflationary models to explain the seeds of 
cosmological perturbations. Measurements on deviations from 
Gaussianity for other primordial bispectral shapes help to shed 
light on more subtle effects about the detailed physics of infia¬ 
tion. 

As in Planck Collaboration XXIV (2014), in the following 
we derive limits on parameters of the models from the NG con¬ 


straints in the following way (unless explicitly stated otherwise): 
we construct a posterior based on the assumption that the sam¬ 
pling distribution is Gaussian (as supported by Gaussian simu¬ 
lations); the likelihood is approximated by the sampling distri¬ 
bution, but centred on the NG estimate (see Eisner & Wandelt 
2009); we employ uniform or Jeffreys’ priors, over intervals of 
the parameters values that are physically meaningful, or as other¬ 
wise stated; and in the cases when two or more parameters are in¬ 
volved, we marginalize the posterior to provide one-dimensional 
constraints on the parameter considered. 


11.1. General single-field models of inflation 

DBI models'. DBI models of infiation (Silverstein & Tong 2004; 
Alishahiha et al. 2004), characterized by a non-standard kinetic 
term of the infiaton field, predict a non-linearity parameter 
= -(35/108)(c“^ -1), where is the sound speed of the in¬ 
fiaton perturbations (Silverstein & Tong 2004; Alishahiha et al. 
2004; Chen et al. 2007b). The corresponding bispectrum shape 
is very close to the equilateral shape. Nonetheless we have con¬ 
strained the exact theoretical (non-separable) shape (see equa¬ 
tion 7 of Planck Collaboration XXIV 2014). The constraint 
we obtain = 2.6 + 61.6 from temperature data = 

15.6 + 37.3 from temperature and polarization) at 68 % CL (with 
ISW-lensing and point sources subtracted, see Table 23) implies 

> 0.069 95 % CL (L-only), (84) 

and 

> 0.087 95 % CL (T+E ). (85) 

In Planck Collaboration XXIV (2014) we constrained the so- 
called infrared (IR) DBI models (Chen 2005b, a), which arise 
in string frameworks. We focused on a minimal setup, consid¬ 
ering a regime where stringy effects are negligible and predic¬ 
tions for primordial perturbations are built within standard field 
theory. In the companion paper Planck Collaboration XX (2016) 
we present an analysis of a more general class of IR DBI models 
which accounts for stringy signatures (see Bean et al. 2008) by 
combining Planck power spectrum and bispectrum constraints. 


Implications for effective field theory of Infiation: In this sub¬ 
section we use the effective field theory approach to infiation in 
order to translate the contraints on f^^^ and f^^^ into limits on 
the parameters of the Lagrangian of general single-field models 
of infiation (of the type P(X, ip) models). In particular we derive 
the most conservative bound on the sound speed of the infiaton 
perturbations for this class of models. 

The effective field theory approach (Cheung et al. 2008; 
Weinberg 2008) provides an efficient way to constrain inflation¬ 
ary perturbations for various classes of models that incorporate 
deviations from the standard single-field slow-roll scenario. In 
this approach the Lagrangian of the system is expanded into 
the (lowest dimension) operators obeying the underlying sym¬ 
metries. We consider general single-field models described by 
the following action 


f 




KHi.2 

-^ Lr 


lidiTiy 


( 86 ) 




-2x. 


+ \mIM\ -c;^) - 


where the curvature perturhation is related to the scalar field n 
as f = -Hn. The infiaton interaction terms 7r(cl,7r)^ and 
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Fig. 29. 68 %, 95 %, and 99.1 % confidence regions in the param¬ 
eter space /nl^°)’ defined by thresholding;^^ as described 

in the text. 

generate two kind of bispectra with amplitudes, respectively, 
/EFTi = -(85/324 )(c- 2 - 1) and = -(10/243)(c-2 - 

1 ) [c 3 + (3/2)Cg j, where M 3 is the amplitude of the operator 

71^ (Senatore et al. 2010, see also Chen et al. 2007b; Chen 
2010b). These two bispectra both peak for equilateral triangles 
in Fourier space. Nevertheless, they are sufficiently different, 
and the total NG signal turns out to be a linear combination of 
the two, leading also to an orthogonal shape. We can put con¬ 
straints on Cs and the dimensionless parameter C 3 (c“^ - 1 ) = 
(Senatore et al. 2010). Notice that DBI inflation¬ 
ary models corresponds to having C 3 = 3(1 - c^)/!, while Cs = 1 
and M 3 = 0 (or C 3 (c“^ - 1 ) = 0 ) represent the non-interacting 
(vanishing NG) case. 

The mean values of the estimators for and are 
expressed in terms of Cs and C 3 by 

1 _ ^2 

/nl”' = —^ [-0.275 - 0.0780c^ - (2/3) X O. 78 OC 3 ] 

1 _ ^2 

^ ^ [ 0.0159 - 0.0167c^ - (2/3) X O.OI 67 C 3 ] . (87) 

Here the coefficients come from the Fisher matrix between the 
equilateral and orthogonal templates and the theoretical bispec¬ 
tra predicted by the two operators and We use a 

statistic given by ;^f^(c 3 , Cs) = p'^(c 3 ,Cs)C"*p(c 3 ,Cs), where 
&(c 3 ,Cs) = /^(c 3 ,Cs) - fp (/^{equilateral, orthogonal}), fp be¬ 
ing the joint estimates of equilateral and orthogonal /nl (see 
Table 11), C the covariance matrix of the joint estimators, and 
/^(c 3 ,Cs) is provided by Eq. (87). As an example in Fig. 29 
we show the 68 %, 95 %, and 99.7 % confidence regions for 
and obtained from the T E constraints, requiring 
< 2.28,5.99, and 11.62 respectively (corresponding to slx^ 
variable with two degrees of freedom). In Fig. 30 we show the 
corresponding confidence regions in the (C 3 , Cs) parameter space. 



Fig. 30. 68 %, 95 %, and 99.7 % confidence regions in the single¬ 
field inflation parameter space (Cs, C 3 ), obtained from Fig. 29 via 
the change of variables in Eq. (87). 


Marginalizing over C 3 we find 

Cs > 0.020 95 % CL (T-only), ( 88 ) 

and 

Cs > 0.024 95 % CL (T-^E ). (89) 

The constraints improve by a few percent in T -only and by up 
to 25 % by including polarization, in comparison with those of 
Planck Collaboration XXIV (2014). 

Galileon models of inflation: Galileon models of in¬ 
flation (Burrage et al. 2011; Kobayashi et al. 2010; 
Mizuno & Koyama 2010; Ohashi & Tsujikawa 2012) are 
well motivated models based on the so called “Galilean sym¬ 
metry” (Nicolis et al. 2009). They are characterized by stability 
properties that are quite well understood (ghost-free, and stable 
against quantum corrections) and can arise naturally within 
fundamental physics models (de Rham & Gabadadze 2010b, a). 
Moreover they are an interesting example of models where 
gravity is modified on large scales and we focus on them 
also as a typical example of a more general class of modified 
gravity theories that are ghost-free (the so called Homdeski 
theories, Horndeski 1974; see also Planck Collaboration XIV 
(2016) where these models are discussed in the context of dark 
energy/modifed gravity scenarios for the late time evolution of 
the universe). The predictions for the primordial perturbations 
are very rich. Bispectra can be generated with the same shapes 
as the “EFTl” and “EFT2” bispectra (see also discussion 
in Creminelli et al. 2011), however the amplitude(s) scale with 
the fluctuation sound speed as differently from the general 
single-field models of inflation considered in the above subsec¬ 
tion. They can be written as (at the lowest-order in slow-roll 
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(90) 


- 2254 + 280- 


(91) 


Here A, Cs, and Cs are dimensionless parameters of the mod¬ 
els. In particular Cs is the sound speed of the Galileon scalar 
field, while Cs is a parameter that appears to break the stan¬ 
dard consistency relation for the tensor-to-scalar perturbation 
ratio (r = 166Cs = being the tensor spectral in- 

dex)^^. Accordingly to Eq. (87) a linear combination of these 
two bispectra generate equilateral and orthogonal bispectra tem¬ 
plates^ ^ From the Planck constraints on and (see 
Table 11), we derive constraints on these model parameters fol¬ 
lowing the procedure described at the beginning of this section. 
We choose log-constant priors in the ranges 10' ^ < A < 10"^, and 
10 “^ < Cs < 10 ^, together with a uniform prior lO""^ < Cs < 1 . 
These priors have been choosen essentially on the basis of per¬ 
turbative regime validity of the theory and to allow for a quite 
wide range of parameter values. In Fig. 31, as an example, prob¬ 
ability contours are shown in the parameter space (Cs,Cs) from 
the T + E constraints, after marginalizing over the parameter A. 
Marginalizing over both A and Cg we find 

cfaiiieon > 0.21 95 % CL (r-only), (92) 

and 

cfalileon > 0.23 95 % CL {T+E). (93) 

Notice that interestingly enough the parameter Cs can be negative 
in principle, corresponding to a blue spectral tilt of infiationary 
gravitational waves (without any kind of instability). We there¬ 
fore also explore this branch, with a log-constant prior (for -Cg), 
- 10 ^ < Cs < - 10 “"^, and the same priors for the other parameters 
as above. Fig. 32 shows the probability contours in the (cs,Cs) 
plane, after marginalizing over the parameter A, for the < 0 
branch. Marginalizing over both A and Cs gives 

^Galileon > q 95 (L-OIlly), (94) 

and 

cfalileon > 0.21 95 % CL (T+E). (95) 

A combined analysis of the Planck bispectrum and power 
spectrum constraints on the Galileon models is presented in the 
companion Planck paper on infiation (Planck Collaboration XX 
2016). 


11.2. Multi-field models 


Curvaton models: The simplest adiabatic curvaton models pre¬ 
dict local NG with an amplitude (Bartolo et al. 2004c,b) 


v^local _ ^ 5 

"4rD 6 3’ 


(96) 


For the explicit expressions of these parameters in terms of the 
coefficients of the Galileon Lagrangian see Planck Collaboration XX 
(2016). 

We note that we are neglecting 0{e\lcl) corrections (where 0{e\) 
means also (9(77s, s,...)); see (Burrage et al. 2011) and (Ribeiro & Seery 
2011). These corrections will have a different shape associated with 
them and they are not necessarily small when compared with some of 
the terms displayed, e.g., the terms (9(1/4). 



Log 10 Cs 

Fig. 31. 68%, 95%, and 99.7% probability contours in the 
Galileon models for Cs and Cs parameters for the Cs > 0 branch 
(tensor spectral index < 0). 



Fig. 32. 68%, 95%, and 99.7% probability contours in the 
Galileon models for c^ and c^ parameters for the Cs < 0 branch 
(blue tensor spectral index > 0). 


for a quadratic potential of the curvaton field (Fyth & Wands 
2002; Fyth et al. 2003; Fyth & Rodriguez 2005; Malik & Fyth 
2006, Sasaki et al. 2006), where rj) — [3/5curvaton/(2/^curvaton 
4pradiation)]D IS the “curvatou decay fraction” at the time of the 
curvaton decay in the sudden decay approximation. Assuming a 
uniform prior, 0 < td < 1, our constraint = 2.5 + 5.7 at 
68 % CF (see Table 11) yields 

td > 0.16 95 % CF (r-only), (97) 

while accounting for temperature and polarization data (/^^^ “ 
0.8 + 5.0 at 68 % CF) gives 

ro > 0.19 95 % CF (T -tE) , (98) 
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improving the previous Planck bound which was previously 
ro > 0.15 (95% CL; Planck Collaboration XXIV 2014). In 
Planck Collaboration XX (2016), assuming there is some relic 
isocurvature fluctuations in the curvaton held, a limit on td is 
derived from the bounds on isocurvature fluctuations. In this re¬ 
stricted case, the limit td > 0.98 (95 % CL) is derived, which is 
consistent with the constraint given here. 

Notice that the above expression of (9^) valid 

under the assumption that there is no signiflcant decay of the 
inflaton into curvaton particles. In general one should account 
for such a possibility. For example, if the classical curvaton held 
survives and starts to dominate, then the curvaton particles pro¬ 
duced during reheating (which have the same equation of state 
as the classical curvaton held) are expected to survive and domi¬ 
nate over other species at the epoch of their decay. The classical 
curvaton held and the curvaton particles decay at the same time, 
inevitably producing adiabatic perturbations (for a detailed dis¬ 
cussion see Linde & Mukhanov 2006). A general formula for 
accounting for the possibility that the inflaton fleld decays 
into curvaton particles, is provided in Sasaki et al. (2006): 


^local 

^NL 


= (1+A2) 


5 

4rD 


5rD 


5 

3 ’ 


(99) 


where = pcurv. particles/Pcurv.fieid Hieasures the ratio of the en¬ 
ergy density of curvaton particles to the energy density of the 
classical curvaton fleld (Linde & Mukhanov 2006; Sasaki et al. 
2006) and Pcurvaton in the expression for td is given by Pcurvaton = 
Pcurv.particies + Pcurv.fieid- Using the Uniform prior 0 < td < 1 and 
0 < Ag < 10^ our measurements of constrain Ag < 8.5 at 
95 % CL (T) and A^ < 6.9 (T-\-E). 


11.3. Non-standard inflation models 

Directional dependence motivated by gauge fields: In Table 25 
we constrained directionally-dependent bispectra (Eq. 23). This 
kind of NG is generated by inflationary models characterized by 
the presence of gauge flelds. An actual realization of this type 
of scenario can be obtained with a coupling between the infla¬ 
ton and the gauge fleld(s), via the kinetic term of the fleld(s), 
i.e., JL = -/^(0)P^/4. In this formula, represents the strength 
of the gauge fleld, while 7(0) is a function of the inflaton fleld 
with an appropriate time dependence (see, e.g., Ratra 1992). In 
this type of scenario, vector flelds can be generated during in¬ 
flation, and this in turn determines the excitation of L = 0 and 
2 modes in the bispectrum, with = XL(\g^\lf).\){Nk^l6Q), 
where = (80/3) and Xi =2 = “(10/6) (Bamaby et al. 

2012b; Bartolo et al. 2013a; Shiraishi et al. 2013a). The param¬ 
eter appearing in the equations above, represents the ampli¬ 
tude of a quadrupolar anisotropy in the power spectrum (see, 
e.g., Ackerman et al. 2007), while N deflnes the number of e- 
folds, before the end of inflation, when the relevant scales exit 
the horizon. It is thus clear that these models predict both a de¬ 
gree of statistical anisotropy in the power spectrum, and a po¬ 
tentially non-negligible bispectrum, as well as a direct relation 
between the two. 

Starting from our SMICA constraints from TiT+E) in 
Table 25, marginalizing over a uniform prior 50 < A < 70, 
and assuming uniform priors on -1 < g^ < 1, we obtain the lim¬ 
its -0.050 < g^ < 0.050 (-0.040 < g^ < 0.040), and -0.31 < 
g^ < 0.31 (-0.29 < g^ < 0.29), from the L = 0, L = 2 modes, 
respectively (95 % CL) (considering g^ as scale-independent). 


We note that these constraints refer to all models in which cur¬ 
vature perturbations are sourced by an fi(fi))E^ term (see ref¬ 
erences in Shiraishi et al. 2013a). The constraints we obtain 
are consistent with the tighter (model-independent) limits ob¬ 
tained in Planck Collaboration XX (2016) for the case of a scale- 
independent g^, from an analysis of quadrupolar anisotropies in 
the CMB power spectrum. 

NGfrom gauge field production during axion inflation: We have 
also constrained the inverse decay NG of Eq. (25) arising typi¬ 
cally in models where the inflaton fleld is a pseudoscalar axion 
that couples to a gauge fleld. Using the modal estimator we get 
the following constraints (removing ISW-lensing bias): 

^mv.dec = 17+43 68 % CL, (100) 

for temperature only; and 

^mv.dec = 23 + 26 68 % CL, (101) 

from temperature-fipolarization (see Table 23). The NG ampli¬ 
tude is given by where 

_ H^/( 27 r\^\) is the power spectrum of vacuum-mode cur¬ 
vature perturbations (i.e., the power spectrum predicted with¬ 
out the coupling to gauge flelds), is the dimensionless 

scalar power spectrum of curvature perturbations (a star denot¬ 
ing evaluation at the pivot scale). The NG parameter is expo¬ 
nentially sensitive to the strength of the coupling between the 
axion and the gauge fleld. From Eq. (101) we limit the strength 
of the coupling to ^ < 2.5 (95 % CL). The details together 

with constraints on the axion decay constant can be found in 

Planck Collaboration XX (2016), where an overview of various 
observational limits on axion (monodromy) models of inflation 
is presented. This limit is in agreement with the one that can be 
derived from tensor non-Gaussianities (see below). 

Tensor NG and pseudoscalars: In inflationary scenarios associ¬ 
ated with a pseudoscalar coupling to a U(l) gauge fleld, the ten¬ 
sor bispectrum generated via the gravitational interaction with 
the gauge fleld is expressed by Eq. (27), and the amplitude 

depends on the following: the coupling strength of the 
pseudoscalar fleld to the gauge fleld {fl)\ a slow-roll parameter 
for the inflaton (6); and the power spectrum of vacuum-mode 
curvature perturbations {P). The expression is ^ 6.4 x 
1011^3^3 (Cook & Sorbo 2013; Shiraishi et al. 2013b). 
The constraints on presented in Table 26 can then be used to 
constrain the model parameters. Clearly there are strong degen¬ 
eracies, but if we marginalize over a uniform prior 1.5 x 10-9 < 
P < 3 X 10-9 and set e = 0.01, then assuming a uniform prior 
0.1 < ^ < 7, from the SMICA (T-only or T-\-E) limit, we obtain 
^<3.3 (95% CL). 

Warm inflation: We update the constraints on warm inflation 
models in the strongly dissipative regime, when dissipative ef¬ 
fects are relevant. In this regime = -151n(l -r /14)-5/2 
(Moss & Xiong 2007) with a large dissipation parameter r^ = 
Tl(3H) (where E is a friction term for the inflaton evolution 
describing the energy transfer from the inflaton fleld to radi¬ 
ation).. The limit from the 2013 Planck release is log^o ^ 
2.6 (95 % CL) (Planck Collaboration XXIV 2014). Assuming 
a constant prior 0 < log^o ^ 4, the new SMICA constraint 
= -23 + 36 at 68 % CL from T (/™^ = -14 + 23 
from T+E), see Table 24, yields a limit on the dissipation pa¬ 
rameter of logio^d ^ 3.3 (logio^d ^ 2.5) at 95% CL, with 
the T+E constraints (in brackets), slightly improving the 2013 
Planck limits. Values of r^ > 2.5 (strongly-dissipative regime) 
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are still allowed; however, the Planck constraint puts the model 
in a regime where there might be an overproduction of graviti- 
nos (see Hall & Peiris 2008 and references therein). Unlike the 
strong dissipative regime, in the intermediate and weak dissi¬ 
pative regimes (ro < 1) the NG level strongly depends on the 
microscopic parameters (T jH and td), giving rise to a new addi¬ 
tional bispectrum shape (for details see Bastero-Gil et al. 2014). 


11.4. Alternatives to inflation 

Ekpyrotic/cyclic models have been proposed as an alternative to 
inflation (for a review, see Lehners 2010). Local NG is gener¬ 
ated from the conversion of “intrinsic” non-Gaussian entropy 
perturbation modes into curvature fluctuations. Models based 
on a conversion taking place during the ekpyrotic phase (the 
so called “ekpyrotic conversion mechanism”) are already ruled 
out (Koyamaetal. 2007; Planck Collaboration XXIV 2014). 
Ekpyrotic models where “kinetic conversion” occurs after the 
ekpyrotic phase predict a local bispectrum with = 

(3/2) a :3 Vi + 5, where the sign depends on the details of the 
conversion process (Lehners & Steinhardt 2008; Lehners 2010; 
Lehners & Steinhardt 2013), where e ^ 50 or greater are typ¬ 
ical values. If we take 6 ^ 100 and a uniform prior on -5 < 
K 3 < 5 the constraints on from T-only (see Table 11), 
yield -0.91 < /C 3 < 0.58 and -0.25 < K 3 < 1.2 at 95% CL, 
for the plus and minus sign in respectively. Lrom the T-\-E 
constraints (Table 11) we obtain -0.94 < /C 3 < 0.38 and 
-0.27 < K3 < 1.0 at 95 % CL, for the plus and minus sign in 
respectively. If we consider 6 ^ 50 we derive the following 
limits: -1.3 < /C 3 < 0.81 and -0.35 < /C 3 < 1.8 at 95 % CL from 
T-only; -1.3 < K 3 < 0.53 and -0.38 < K 3 < 1.5 at 95 % CL from 
T+E. Another variant of the ekpyrotic models has been investi¬ 
gated in (Qiu et al. 2013; Li 2013; Lertig et al. 2014), where the 
intrinsic NG is zero and NG is generated only by non-linearities 
in the conversion mechanism, reaching a value of ^ +5. 


11.5. Inflationary interpretation of CMS trispectrum results 


We briefly interpret the trispectrum constraints in an inflation¬ 
ary context. Lirst, consider the case of single held inflation. The 
action for the Goldstone boson tt is highly constrained by resid¬ 
ual diffeomorphism invariance (see, e.g.. Smith et al. (2015)). To 
lowest order in the derivative expansion, the most general action 
is: 






+ 2 M^ 


-Vi? - + {dn7T?(dy7T? 


—^[^kUnk\d^nf + ---] 

+ ^[\6kU32k\d,nf+ 


( 102 ) 


where the parameter M 4 is related to the trispectrum by 


288 ^ ® 


(103) 


The constraint in Eq. (70) translates to the following param¬ 
eter constraint in single held inflation: 


- 9.70 X 10*^ < 




H^cl 


< 8.59 X 10'"^ 


(95% CL). (104) 



Fig. 33. 68 % and 95 % confidence regions in the ) 

plane, with the Lorentz invariant model in Eq. (107) shown as 
the dashed line. 


This constraint is a factor 1.8 better than WMAP. 

Turning now to multifleld inflation, we consider an action of 
the form 


So- = J"d^x 


1.9 1 4 

-(dnCr)^ -\ - 

/4 9 ^4 




(105) 


where cr is a spectator held that acquires quantum fluctuations 
with power spectrum Pcr(k) = I( 2 k^) and converts to adiabatic 
curvature via f = (2A^y^^H~^cr. The trispectrum in this model 
is 



^NL 


A, = 




A, = 


25 

325 

2575 774 
20736X4’ 


(106) 


so we can constrain its parameters by thresholding the;^^ deflned 
in Eq. (75). Lor example, if we consider the Lorentz invariant 
model 



ziducrf 




(107) 


so that the parameters A/ of the more general action in Eq. (105) 
are given by Af = - 2 A 2 = A^ = A^, then by thresholding at 
= 4 (as appropriate for one degree of freedom), we obtain 
the following constraint on the parameter A: 

- 0.26 < ^ < 0.20 (95 % CL). (108) 

Constraints in other parameter spaces can also be obtained by 
thresholding the;^^ deflned in Eq. (75). Lor example we show 

68 % and 95 % confldence regions in the (^^L’ )-plane, ob¬ 
tained by thresholding at;^^ = 2.28 and;^^ = 5.99 as appropriate 
for a;^^ random variable with two degrees of freedom. 
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DBI Trispectrum: The trispectrum constraint on the shape in 
Eq. (70) can also be used to obtain a lower bound on the DBI 
model sound speed. This is because, in the small sound speed 
limit (Chen et al. 2009; Arroja et al. 2009), the dominant con¬ 
tribution to the contact interaction trispectrum (Huang & Shiu 
2006) has this shape. The corresponding non-linearity param¬ 
eter is = -25/(768 c^). We follow the same procedure as 
described at the beginning of this section 11 and, assuming a 
uniform prior in the range 0 < < 1/5, we can derive a con¬ 

straint on Cs as 

c™’> 0.021, 95% CL. (109) 

This constraint is consistent with the ones derived from the bis¬ 
pectrum measurements (see Eqs. 84 and 85)) and it is only a 
factor of about three worse. Notice, however, that in this case 
we are ignoring the scalar exchange contribution, which is of the 
same order in Cg. 

Curvaton trispectrum: Eor the simplest curvaton scenario, 
the trispectrum non-linearity parameter prediction is 

(Sasaki et al. 2006) 

^ 2 

Eollowing the procedure described at the beginning of Sect. 11, 
we use the observational constraint obtained in Sect. 9 (Eq. 70), 
and the same prior (0 < td < 1) as in 11.2, to obtain a lower 
bound on the curvaton decay fraction as 

rD> 0.05 95% CL. (Ill) 

This limit is consistent with the previous ones derived using the 
bispectrum measurements and it is a factor of about 3 to 4 worse. 

12. Conclusions 

In this paper we have presented the constraints on pri¬ 
mordial NG using the full Planck mission data. The re¬ 
sults have improved compared to the Planck 2013 re¬ 
lease (Planck Collaboration XXIV 2014) as a consequence of 
including data from the full mission and taking advantage of 
Planck's polarization capability — the first time that maps of 
the CMB polarization anisotropies have been used to constrain 
primordial NG. 

Using temperature data alone, the constraints on the local, 
equilateral, and orthogonal bispectrum templates are “ 

2.5 ± 5.7, = -16 ± 70, and = -34 ± 33. Moving from 

the nominal Planck 2013 data to the full mission data yielded 
modest improvements of up to 15 % (for the orthogonal shape). 
After the inclusion of full mission polarization data, our final 
constraints become = 0.8 + 5.0, = -4 + 43, and 

±21, which represents a substantial step forward 
relative to Planck 2013, with error bars shrinking by 14% (lo¬ 
cal), 43 % (equilateral), and 46 % (orthogonal). These improved 
limits on the standard shapes enhance our understanding of dif¬ 
ferent inflationary models that can potentially lead to subtle ef¬ 
fects beyond the simplest models of inflation. 

The reason that the polarization data provide such comple¬ 
mentary constraints on primordial curvature perturbations is due 
to the phase shift of the CMB polarization transfer functions 
compared to the temperature transfer functions. So, despite the 
relatively much lower signal-to-noise in the polarization maps, 
their inclusion leads to appreciable improvements on limits on 


NG parameters. Nevertheless, the full characterization of the 
noise properties in the polarized maps is still ongoing. In spite of 
the extensive testing and cross-checks validating the combined 
temperature and polarization results, we therefore conservatively 
recommend that all results that include polarization information, 
not just the polarization-only results, be taken as preliminary at 
this stage. 

The complementary nature of the polarization information 
also represents an important cross-check on the analysis. The 
Planck results based on polarization alone are statistically con¬ 
sistent with the results based on temperature alone, with a pre¬ 
cision comparable to that achievable in an optimal analysis 
of the WMAP 3-year temperature maps (Spergel et al. 2007; 
Creminelli et al. 2007; Yadav & Wandelt 2008). 

Our analysis was subject to an extensive validation exercise. 
In addition to extensive simulation tests, including for the first 
time a detailed test of the impact of time-domain de-glitching, 
our results are supported by tests for robustness under change of 
estimator implementations (KSW, binned bispectrum, and two 
modal estimators), and variations in sky coverage as well as up¬ 
per and lower multipole cutoffs. We also tested for possible di¬ 
rectional dependence using a needlet estimator. These tests form 
the basis of our selection of SMICA as the main foreground clean¬ 
ing method for our headline results. 

The Planck 2015 analysis presented here provides con¬ 
straints on a greatly extended range of template families. These 
extensions include a tenfold increase in the range of frequen¬ 
cies covered in feature models, giving rise to linearly oscillat¬ 
ing bispectra, generalized shapes for oscillating models includ¬ 
ing for logarithmic oscillations, tests for deviations from the 
Bunch-Davies vacuum, models of equilateral type in the context 
of the effective field theory of inflation, and direction-dependent 
primordial NG. Beyond purely scalar mode templates we also 
tested for parity-violating tensor NG. 

Using the full mission data with polarization, we have in¬ 
vestigated the hints of NG reported in the Planck 2013 analy¬ 
sis of oscillatory features. While no individual feature or res¬ 
onance model rises above our detection threshold of 3 cr (after 
inclusion of the look-elsewhere effect), the results of integrated 
(multi-peak) statistical tests indicate that continued investigation 
of oscillatory and non-scaling models is warranted. 

In addition to searches for specific NG templates, we present 
model-independent reconstructions of the temperature and po¬ 
larization bispectra using the modal and binned bispectrum ap¬ 
proaches. These full mission reconstructions can achieve twice 
the resolution of the Planck 2013 results, demonstrating excel¬ 
lent consistency in temperature, and good agreement with the 
WMAP9 reconstruction in regions where this earlier data set is 
signal-dominated. 

The inclusion of polarization information leads to signifi¬ 
cantly improved constraints on NG in primordial isocurvature 
perturbations, providing complementary information to 2-point 
function results for models where the NG in isocurvature compo¬ 
nents is more easily detectable than its contribution to the power 
spectrum. 

A significant addition to this year’s analysis is the inclu¬ 
sion of detailed trispectrum results due to cubic NG. The lo¬ 
cal trispectrum is constrained by Planck temperature data to be 
“ (-9.0 ± 7.7) X lO'^ and the other two shapes were also 
found to be consistent with Gaussianity. Both 3-point and 4- 
point constraints are consistent with the improved (though still 
suboptimal) constraints from Minkowski functionals, a very dif¬ 
ferent estimation framework. This concordance adds confidence 
in our results. 
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We have discussed the implications of our results on the 
physics of the early Universe, showing that the ^-point func¬ 
tions for ^ > 2 provide a significant window onto the primor¬ 
dial Universe beyond the power spectrum, constraining general- 
single field, multifield, and non-standard infiation models, as 
well as alternatives to infiation. Using bispectrum and trispec¬ 
trum limits we updated results on the parameter space of the in¬ 
fiationary models (and alternatives) already tested in 2013, and 
constrained the parameter space of other well-motivated infla¬ 
tionary models (e.g., Galileon-like models of infiation, and mod¬ 
els where axion/pseudoscalar fields are present during infiation). 

The global picture that emerges is one of consistency with 
the premises of the ACDM cosmology, namely that the structure 
we observe today is the consequence of passive evolution of adi¬ 
abatic, Gaussian, primordial seed perturbations. Nevertheless, 
NG at some level is expected in all infiationary models, and 
hence we should strive to find means to reduce errors on /nl 
still further. 
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Appendix A: Derivation of an estimator for cl 

As parameterized by Eq. (23), we express a primordial bispec- 
tmm of direction-dependence: 


B^(ki,k2,k3) = 


^ Cl [Piih • k 2 )P^(ki)P<^(k 2 ) + 2perm.] 
L>1 


(A.l) 


where Piiki • ^ 2 ) is a Legendre Polynomial of order L. It can 
be shown that such a primordial bispectrum leads to a CMB bis¬ 
pectrum 
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where p denotes either temperature or £’-mode polarization, 
is the reduced CMB bispectrum, and the term with 


large parentheses denotes the Wigner-3j symbol. The reduced 

iIclXpi, 


CMB bispectrum is given by (Shiraishi et al. 2013a) 


i1cl),P\P2P2 


An H r^dr 

2L+ 1 Jo 
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(A.3) 


where denotes the Kronecker delta function, {...} the Wigner 
6j symbol, “perm.” means permutations, W£ is a beam window 
function, and 

dink g^^(k)j((kr), (A.4) 

= “/ dC\nkPg,{k)g’’^{k)jt'{kr). (A.5) 

Here, (k) is the radiation transfer function for temperature or 
£’-mode polarization, and j^{x) is a spherical Bessel function. 
The h symbol is 




(2Ai -f \){2h + 1)(2^3 + 1) 
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By maximizing the likelihood with respect to cl, we obtain the 
KSW estimator for 

- _ ^ ^rnim2m3AcL),P\P2P3 

^ ~ 6 N ^ ^ 

PiQi iinii 

-3((C-‘«C(C-'<„,)mc (C-'<„3l, (A.7) 

where MC denotes that the average is over Monte Carlo simu¬ 
lations and TVc^ is a normalization constant. The normalization 
constant Nc^ is given by 

(A.8) 

PiQi 

^ y.CL),P\P2P3 ^^-\^V2q2 r(y-l)^3^3 ^(Cl),^1^2^3 

Using Eq. (A.2), we find that 
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Since b^!^Xr) ^ (-1)"^ ’ ^Luin, r) is not a real func¬ 

tion, but a complex function. We estimate Bimin, r) efficiently 
by computing the following with HEALPix (Gorski et al. 2005): 


Re[BtM(n,?-)] = (r) Y(„{ny, 

im 

£m 

Here 


r.. b^’^{r) + {-\T[b^’^ {r)Y 
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Appendix B: Definition of Minkowski Functionais 
and theoreticai expectations 

For a field f(x) of zero average and variance cr^ defined on the 
two-dimensional sphere S^, an overdense excursion set is de¬ 
fined as 

Z = {xeS^\f(x) > vcro}. (B.l) 

The boundary of the excursion is 

(9Z = {x G S^l f(x) = vo-o). (B.2) 

Then the three Minkowski functionals on the sphere are 


where 

b(,Zir) = ht'CL 
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withbeing defined in Eq. (A.5). In the derivation above, 
we used the identities 
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and 

6(ki + k2 + ki) = S. f drdn^ J^/'jy(kir)Ye,„,(ki)ylm,(k) 

^ Ami 

X Z Z i'^h(k3r)Ye,„,,(ki) 

A m2 A m3 

Applying Eq. (A.9) to Eq. (A.7), we find 


where df2 and dl are respectively elements of solid angles (sur¬ 
face) and of angle (distance), k is the geodesic curvature. Note 
that the Genus can be also expressed as the number of compo¬ 
nents^^ in the excursion minus the number of holes in the excur¬ 
sion. 

The fourth functional y3(y), is defined, for y > 0, as the 
number of components in the excursion. Symmetrically, for y < 
^^0, it is the number of underdense components (or the number of 
components in the excursion {x e S^| f(x) < vctq}). 

In the Gaussian limit, the functionals can be expressed the 
following way (see, e.g., Matsubara 2010; Vanmarcke 1983): 
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A component is a connected subset of the excursion. 
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The amplitude depends only on the shape of the power spec¬ 
trum Q: 




^3 
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Cl>2 


\2 


2/ cri 

V20-0 



k<2 
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(B.ll) 


where cok = 7r^/^/r(^/2-r 1), which gives ojq = 1, tui = 2, (jl >2 = n 
and (To and (T\ are respectively the rms of the field and its first 
derivatives. 
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